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ABSTRACT 


Effective analysis tools have been developed for predicting the nonlinear 
rotordynamic behavior of the SSME turbopumps under steady and transient operating 
conditions. Using these methods, preliminary parametric studies have been conducted on 
both generic and actual HPOTP (high pressure oxygen turbopumps) models. In particular, 
a novel modified harmonic balance/alternating Fourier transform (HB/AFT) method was 
developed and used to conduct a preliminary study of the effects of fluid, bearing and seal 
forces on the unbalanced response of a Multi-disk rotor in presence of bearing clearances. 
A computer program was developed and made available to NASA, Marshall. The method 
makes it possible to determine periodic, sub-, super- synchronous and chaotic responses of 
a rotor system. The method also yields information about the stability of the obtained 
response, thus allowing bifurcation analyses. This provides a more effective capability for 
predicting the response under transient conditions by searching in proximity of resonance 
peaks. Preliminary results were also obtained for the nonlinear transient response of an 
actual HPOTP model using an efficient, newly developed numerical method based on 
convolution integration. A computer program was developed and made available to NASA 
Marshal Flight Center. Currently, the HB/AFT is being extended for determining the 
aperiodic response of nonlinear systems. Initial results shows the method to be promising. 
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I. INTRODUCTION 


Background : 

Modern mechanical systems are being recently designed for higher performance, 
reliability and smooth operation within compact configuration. These requirements often 
cause significant nonlinear effects which could not be predicted with linear models. 
Therefore, a more complete picture of their nonlinear dynamic characteristics is required 
to enhance their efficient design, refinement, monitoring or maintenance. 

Modern complex rotating machinery, such as the turbopumps of the space shuttle 
main engines (SSME), contain various sources of strong nonlinearities. These include 
clearances and nonlinearity of rolling elements, rubbing in splines and in built-up rotor 
segments, rubbing at seals and rotor blades, viscous damping and various fluid effects. 
Observed nonlinear behavior of actual rotor systems include jump discontinuities [1], large 
subsynchronous motion, [2] - [4], quasi-periodic and possible chaotic motion [5]. As stated 
by Nataraj and Nelson [6], the future developments in modern machines heavily depends 
on the ability to identify, understand mathematically and analyze systems involving nonlinear 
components. This is particularly the case for the proper development, monitoring and 
analysis of the SSME turbopumps. 

Quite often, it is essential to determine the steady state periodic response of rotor 
systems in the form of self excited limit cycles or forced motion due to rotating imbalance. 
Accurate prediction of the nonlinear periodic responses and their stability plays a central 
role in developing a complete picture of the dynamic behavior of nonlinear rotor systems 
as function of their parameters. 
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Several methods have recently been advanced for determining the periodic response 
of low order nonlinear rotor systems, [7] - [10]. For application to large, multi-disk rotor 
systems, Nataraj and Nelson [6] developed a periodic solution method based on a 
collocation approach for the response of the rotor. They utilized a subsystem approach to 
reduce the size of the resulting system of algebraic equations. Ehrich [11] recently analyzed 
high order subharmonic response and chaos [12] using numerical integration for Jeffcott 
rotor with a bearing clearance. 

Few analysts have addressed the stability of periodic or subharmonic responses of 
nonlinear rotor systems despite its considerable significance in predicting the response of 
modern, high performance systems. Most of the stability or bifurcation analyses were 
concentrated on one dimensional (rectilinear) problems where solution forms are assumed 
a priori (Shaw and Holmes [13], and Natsiavaas [14]). These approaches are very difficult 
or not feasible to extend to the two dimensional nonlinear rotor problems. 

Complete characterization of the dynamic behavior must include determining the 
steady state responses and their bifurcation as function of the system parameters. In 
addition, the transient response of the nonlinear system has to be determined as part of any 
complete dynamic analysis of the system. 

Objectives and Outline of Study 

The main objective of this study was to develop reliable and efficient analytical- 
computational methods of the nonlinear dynamic analysis of large, rotor-housing systems 
such as the turbopumps of the space shuttle main engines (SSME), and some aspects were 
then examined to be of the nonlinear behavior of a general multi-disk rotor-housing system. 
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In the present study, HB (Harmonic Balance)/AFT (Alternating Frequency-Time) 
method has been developed. Using the method, a study is made of the dynamics and 
stability of simplified models of the HPOTP (High Pressure Oxygen Turbo Pump) of the 
SSME (Space Shuttle Main Engine) including clearances between the bearings’ outer races 
and their supports. The method employs an explicit Jacobian form in an iterative procedure 
which ensures convergence at all parameter values. A dynamic reduction technique [15] is 
used with the HBM to reduce the system nonlinear differential equations to linear algebraic 

equations involving only the nonlinear coordinates. 

A recently developed convolution approach by the present author for the non-linear 
transient analysis was carefully tested against direct integration techniques and proved more 
robust and efficient. Using the HB/AFT method, the response resonances of a multiple- 
disk rotor as function of the rotor spinning speed (critical speeds) can be located. The 
convolution approach can then be used effectively to determine the transient response in 
passing through these critical speeds. 
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II. ANALYSIS METHODS AND RESULTS 


STEADY STATE RESPONSE AND STABILITY 
(i) One and Two Dimensional Systems: 

The dynamic behavior of strongly nonlinear mechanical systems with piecewise-linear 
or piecewise-smooth nonlinearities is studied using a newly developed HB (Harmonic 
Balance) with an AFT (Alternating Frequency Time) robust and efficient algorithm. By 
employing the harmonic balance approach, the nonlinear differential equations are 
transformed to nonlinear algebraic equations. Two iterative techniques are available for 
solving the nonlinear algebraic equations. These are the Newton-Raphson method and a 
certain version of the Quasi-Newton Algorithm. The Newton-Raphson algorithm requires 
a laborious and complex calculation of a Jacobian matrix and has a narrow range of initial 
guesses for achieving convergence. However, the algorithm has superior convergence speed 
as compared to the corresponding Quasi-Newton method. 

The development of the method and its application to an oscillator interacting 
through a gap with a flexible stop (see Figure 1) can be found in a paper [16], accepted for 
publication in the Journal of Applied Mechanics of the American Society of Mechanical 
Engineers (ASME). The method was also applied to a modified Jeffcott rotor model (see 
Figure 2) supported on bearings, with clearances (see appendix [A] which was also published 

as a paper in "Nonlinear Dynamics", [17]). 

A bifurcation analysis method, based on Flouqet’s theory, was also developed for 
determining the stability of the obtained periodic solutions. For the stability analysis, 
Poincare mapping is utilized to obtain the fixed points corresponding to the periodic 
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Figure 1-a Oscillator with a gap 


s(*) 



Figure 1-b Piecewise-1 inear restoring force of 
oscillator 
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solutions (or limit cycles). Small perturbation around the periodic solutions (fixed points) 
is performed in order to analyze their stability. A first order Jacobian around a fixed point 
is calculated using numerical integration to obtain the associated monodromy matrix. The 
eigenvalues of the monodromy matrix are analyzed to determine the bifurcation type (cyclic 

fold, secondary Hopf or period multiplying). 

The bifurcation analysis allows determining ranges of parameters at which the 
response of a given rotor system would become subsynchronous or chaotic. 

(ii) Multi-Disk Rotor Systems 

The HB/AFT method was also applied to a two-disk rotor system (see Figure 3) 
containing a bearing clearance [18,19]. The method generalizes the author s early work 
[20,21]. Results were obtained for the synchronous, sub-synchronous and chaotic response 
of the system. The bifurcation analysis developed for this case was also used for predicting 
the onset of qualitative changes in the dynamic behavior of the system. The extended 
HB/AFT and bifurcation analysis extended for application to the two-disk rotor system 
considered, as well as the sample results obtained, are included in Appendix B. The work 
in this appendix was also accepted for publication in "Nonlinear Dynamics Journal [18]. 
A computer program for this case was developed and made available to NASA, Marshall 
which, if needed, can be readily modified for modeling the turbopumps of the SSME. 

NONLINEAR TRANSIENT ANALYSIS 

A convolution approach first reported by Noah [22] and Noah, et. al. [23] was further 
refined and applied to the transient analysis of the generic model shown in Figure 4. A 
general purpose computer program is written based on this approach, along with a user 


7 



Figure 3 . Multi-disk rotor model with a bearing clearance. 
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ROTOR: Shaft diameter: OD = 3.0 inches, ID = 0.0 inch 

Material: Steel, E = 3.0 x 10* psi 
Joint length = 3.0 inches 
Rotor length = 27.0 inches 

HOUSING: Housing/rotor weight ratio = 6/1 

(EI)r = 1.1928 x 10* lb m s 
(EI)h = 2-4 x 10* lb in 5 

Fi(t) = mj e<^> 2 cos4>t hsz, — 4.0 x 10* lb/ in 

Fi(t) = m.\e4> 7 sin4>t Esy, = 5.0 x 10* lb /in 

Fi(t) — mj e<£ 3 cos«f>t Eszj = 5.0 x 10* lb /in 

F 4 (t) = mj e<ji 3 sintpt A’sVj = 1*5 x 10* lb /in 

F s (t ) = mje^cos^t K a = 5.0 x 10 s lb /in 
Ft(t) = mje^ 1 sin</>i Cq = 0.0 


Figure 4. The generic model of the SSME turbopurtp 
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manual and was provided to NASA, Marshall. 

The convolution approach can be applied to a locally nonlinear general rotor housing 
system with rotor imbalance during start-up or shut-down. In the present work, eigen- 
coordinates are used to represent both housing and rotor. The local nonlinearities were 
taken as bearing deadband clearances at the rolling element bearings which support the 
rotor in its housing. The integral formulation of the rotor motion is represented by its 
transition matrix and that of the housing by a convolution integral (based on the housing s 
impulse response). 

The convoluted impulse response can only be applied to a system of uncoupled 
equations while the transition matrix formulation, in addition, can be applied to coupled 
equations. The transition matrix can therefore be applied to coupled dynamical systems 
represented by their physical coordinates or, in case of rotors, coupled by the gyroscopic 

terms in otherwise decoupled modal representation. 

Sample of the various tests conducted on the convolution method, to test its accuracy 
and efficiency, is presented in Figures 5 and 6 for the generic model of Figure 4 under 
transient loading. It can be seen that the method is more robust and efficient than direct 
numerical integration. 

More complete presentation of the method and the transient response is included in 
Appendix C . A paper [24] was also published based on this appendix in the ASME Journal 
of Applied Mechanics. 


10 


500 -i 


LINEAR 


♦ 


Q 

Z 

o 

o 

LJ 

{/) 


U 

2 


z> 

CL 

o 


400 A 


DUHAMEL INTEGRAL 
transition matrix 

- * - RUNGE KUTTA 

— 0- - hybrid 


300 A 


200 A 


loon 


o+- 

0.0 


i ♦ 

' / 

/ 



TIME STEP SIZE (X 1000 SECOND) 


1. 
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III. CONCLUSIONS AND RECOMMENDATIONS 


SUMMARY AND CONCLUSIONS 
(i) Steady State Analysis; 

A robust iterative numerical procedure based on the HBM/AFT method has been 
developed for obtaining the periodic response of a large rotor/housing system containing 
bearing clearances. Modern bifurcation theory is utilized to characterize the dynamic 
behavior of the system. A bifurcation analysis method is developed which provides 
boundaries of parameter regions at which rotor whirling pattern changes its shape rapidly, 
resulting in the occurrence of subharmonic, aperiodic or possible chaotic motion. 

Results on the effects of parameters on a SDOF (Single Degree of Freedom) system 
with piecewise-linear response show that, for some combinations of these parameters, the 
system response exhibits both period doubling and saddle-node bifurcations. Chaotic motion 
was also observed for finite stiffness ratios. The stability analysis, along with the harmonic 
balance-based method provide a very powerful tool for better understanding of the behavior 
of systems with clearances. 

The main results obtained using the HBM/AFT approach with a nonlinear Jeffcott 
rotor 1 with a bearing clearance can be summarized as follows [17]: 

1. Increasing the dimensionless support stiffness ratio, «, causes flip bifurcation 
to occur which produces period doubled whirling motion (subsynchronous 

motion). 


l A Computer program based on this work was submitted to NASA, Marshall. 
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2. For the same a and a nondimensional frequency rut io. d> ■ an increase in the 
damping ratio, c, leads to elimination of the subharmonic motion. 

3. With higher values of a, chaotic whirling motion is feasible. 

4. Limited results obtained concerning the effect of Coulomb friction indicate 
that the coefficient of friction, n, has little effect on the subharmonic 
response. However, higher fi could eliminate the subharmonics near flip 
(period doubling) bifurcation boundaries. 

5. Increasing the nondimensional cross coupling stiffness ratio y leads to a Hopf 
bifurcation which could result in aperiod i c whirling. 

A further developed HBM, using a DFT (Discrete Fourier Transform)/IDFT 
(Inverse Discrete Fourier Transform), is employed to obtain the steady-state periodic 
response for MDOF (multi-degree of freedom) rotor systems with bearing clearances 
(piecewise-linear type nonlinearity). 2 A dynamic reduction (impedance) technique [15] is 
utilized to reduce the system to only those of the nonlinear coordinates. The stability 
analysis is performed via perturbation of the obtained periodic solutions. The reduced and 
approximated system parameters (mass, damping, and stiffness) are calculated from the 
determined harmonic coefficients. 

A simple MDOF rotor system with a bearing clearance is used for illustration of the 
method. The results obtained show that: i) the HB/AFT method, as developed here, is 
robust and efficient; ii) the method leads to accurate bifurcation boundaries for nonlinear 
MDOF rotor systems and, furthermore; iii) the method can in general be applied to MDOF 
rotor systems with piecewise-smooth or polynomial type nonlinearities at the bearing 

2 A computer program based on this work was submitted to NASA, Marshall. 
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supports. 


In summary, the major advantages of the HB/AFT approach are: 

a) it can provide steady-state periodic solution as well as steady-state quasi- 
period solutions using modified DFT/IDFT algorithms. 

b) its formulation is neither problem nor response pattern dependent, except for 
selecting the least appropriate numbers of harmonics for the Fourier 

expansions. 

c) it can drastically reduce the computational time while providing high 
computational accuracy. Especially for nonlinear MDOF systems with lower 
damping, this approach is much more powerful in delineating the steady-state 

solutions efficiently. 

d) it enables perturbation of the determined periodic solution so that the 
resulting linear ordinary differential equation with periodic coefficient would 
yield information about the stability employing Floquet theory. Approximate, 
but accurate, stability information can then be obtained. 

e) using modern dynamical theory, detailed bifurcation boundaries and their type 
(such as flip, Hopf and fold) can be easily calculated from the perturbed 

periodic solution as function of the system parameters. 

f) by observing the unstable solutions, bifurcation boundaries can be easily 

obtained as function of system parameters. 

(ii) Transient Analysis: 

The hybrid convolution approach, further developed in this study, is shown to provide 
an efficient and accurate closed form integral formulation for determining the transient 
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response of linear systems coupled through local nonlinearities associated with friction and 
clearances. A typical application in which the present method proved quite effective is the 
determination of the transient response of a generic model of the high pressure oxygen 
turbopump (HPOTP) of a space shuttle main engine (SSME) in presence of bearing 

clearances (see Appendix C ). 3 

The use of the transition matrix is successful in allowing the representation of rotor- 
system involving skew-symmetric matrices of gyroscopic loads or other nonconservative 
systems with general velocity dependent matrices. A convolution integral would represent 
quite effectively other systems with normal modes, such as the housing of the HPOTP or 
other non-rotating, proportionally damped structures. An application is also made to a 
generic model of the high pressure oxygen turbopump (HPOTP) of a space shuttle mam 
engine (SSME) in the presence of bearing clearances, constituting the local nonlinearities. 

Two iterative techniques, the Jacobi method and the Gauss-Seidel method, were 
studied. Both were able to correct the predicted coupling forces and converge to the correct 
force magnitudes with a desired accuracy. The Gauss-Seidel method is more efficient in 
CPU time especially for a system with large time step and large external forces. For a 
system with very small external forces, the advantage will not be significant. However, the 
solution has to be formulated differently so as to accommodate a given type of nonlinear 
component. If the convolution method is applied to a system with a specific nonlinear 
component only, Gauss-Seidel scheme could be a better choice. On the other hand, the 
Tamhi scheme is more flexible in its application. Once the equations of a rotor and its 
housing are derived, they are ready for use with almost all other types of coupling 


3 A computer program 


based on this work was submitted to NASA, Marshall. 


components. 

A generic model of the HPOTP of the SSME was used to test the transient method 
and conduct parametric studies, typical of nonlinear rotor/housing systems with bearing 
clearances. 

As typical of the effectiveness of the convolution method, firstly, the accuracy and 
CPU time are studied using both the hybrid convolution approach and Runge Kutta 4th 
order method. The results show for the given triangular load used in [24] and included in 
Appendix C that: 

1. For a fixed tolerance (1 x 10*), the hybrid convolution method is faster and 
more accurate. The CPU time for the Runge Kutta is 1.42 times of the 
hybrid method. For more meaningful comparison, the accuracies of the 
Runge Kutta and convolution methods are made closer by reducing the 
allowable tolerance for the Runge Kutta to 1 x 10 14 . The CPU time for the 
Runge Kutta increases quickly from 1.42 times of the hybrid method s CPU 

to 4.23 times. 

2. The hybrid convolution method is also more robust than the Runge Kutta 
method. The Runge Kutta algorithm failed to converge for time increments 
greater than 2 x KT 5 seconds. However, the hybrid method will diverge when 
the time step size is larger than 3.3 x 10' 5 seconds. 

Secondly, the convolution method is shown to have the potential as a useful tool for 
determining transient response. For the generic model of the HPOTP, three nonstationary 
cases have been studied [25]. One is a linear model with no gap. In this case, a closed form 
solution can be used directly. The other two are nonlinear cases with small and large gaps 
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(gap size 5 x 10 J and 1 x 10' 3 inches respectively). The gaps are considered as the small and 

large bearing clearances. The study shows: 

1. The reaction force in the deceleration period is much higher than that in the 

acceleration period. One of the reasons is that the response in the 
deceleration period has more time to build up since the absolute value of the 
deceleration is smaller than that of the acceleration. 

2. The small gap reduces the approximate first linear critical speed by 9%. 

However, the large gap reduces the speed by 18%. 

3. The small gap reduces the amplitude of the peak bearing forces at the first 
critical speed of the linear case by 45.2% in bearing 1 and 30.6% in bearing 
2. The large gap reduces the amplitude of the peak bearing forces at the first 
critical speed of the linear case by 49.8% in bearing and 34.9% in bearing 2. 
The larger gap causes more reduction in peak forces. 

4. In contrast to the first critical speed, the second critical speed and its 
corresponding peak forces in the two bearings in the deceleration period of 
these three cases are not significantly influenced by the existence or the size 
of the gap. The average reductions are 0.9% for critical speeds and 2.4% for 

the corresponding forces. 

5. The peak force at the second critical speed is higher than that at the first 
critical speed. The ratio is 1.78 and 3.94 times in bearings 1 and 2, 
respectively, for the linear case. For the nonlinear case, the ratios are much 
higher. For the small gap, the ratios are 3.4 and 5.8 in bearing 1 and 2, 


17 


respectively. For the large gap, the ratios are 3.7 and 6.2 in bearing 1 and 2, 

respectively. 

The convolution formulation allows accommodating with ease changes in the 
nonlinear or linear coupling parameters among the various linear subsystems involved. 
Besides bearing clearances, other cases of local nonlinearities involving dry friction and 
impacts were also studied in [25]. 
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RECOMMENDATIONS 


The newly devised HB/AFT method proved to be effective in obtaining the steady 
state solutions for multi-disk rotor systems. Although the impedance (dynamic 
condensation) method [22,15] is utilized to reduce the systems’ equations, the reduction was 
only possible in physical coordinate systems. In order to adopt the method to generalized 
coordinate systems, modal representation utilizing a drastically truncated set of modes will 
be necessary. Further development of the HB/AFT method should be made in order to 
allow for better reduction techniques in modal coordinates. In addition, the so called 
internal resonances (in which nonrational relation between the various modes exists) need 
to be investigated as related to the turbopumps of the SSME. 

Further study is needed to complete the development of a preliminary method 
obtained in this study concerning quasi-periodic response, and to generalize the HB/AFT 
method to conduct parametric studies on the SSME turbopumps under various operating 
conditions. 

Modern studies have revealed the significance of predicting what is labeled "crisis". 
This is a generalization of the jump behavior in nonlinear systems. A preliminary study was 
made and limited results were obtained during the course of the present study. It is 
recommended that further study be pursued in which the "crisis behavior" in piecewise- 
smooth systems (specially in applications to rotordynamics) is peculiar. In this connection, 
the HB/AFT method could prove to be very effective. 

For the transient analysis, further work should be made on the analysis and behavior 
of rotor systems with bearing clearances and rubs. This could include the following: 

a. Use approximate methods to replace local nonlinearities with linear 
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components and then use closed form convolution representation for the 
solutions. 

b. Use other possible iterative techniques, possibly incorporating the Gauss- 
Seidel technique. 

c. Develop other means of increasing the efficiency of the method, including the 
use of predictor-corrector and other algorithms. 

d. Adapt the HB/AFT method to determine domains of attractions of the 

nonlinear rotor system. 

e. Extend the method for application to other local nonlinearities, e.g. Coulomb 
friction rubs at seals, turbine blades, rotor shrink fits and to impacts due to 
intermittent contacts. 

f. Conduct parametric studies of the turbopumps of the SSME using the current 
and modified versions of the computer programs provided to NASA, Marshall, 
developed in this study to examine the significance of the various nonlinear 
phenomena. 

g. Explore the HB /AFT method role in experimental verification, identification, 
and monitoring of the SSME systems. 

Finally, direct comparison should be made of the techniques reported here to those 
currently used throughout industry. It is anticipated that the techniques described here will 
yield more information about the dynamic behavior of nonlinear rotor systems in an efficient 

and systematic manner. 
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"Bifurcation Analysis for a Modified Jeffcott Rotor with Bearing Clearances," Kim, Y.B. and 
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Bifurcation Analysis for a Modified Jeffcott Rotor with Bearing 
Clearances 


Y. B. KIM and S. T. NOAH 

Mechanical Engineering Department , Texas A&M University, College Station. TX 77843 -3123, U.S.A 


Abstract A HB (Harmonic Balance) /AFT (Alternating Frequency /Time) technique is developed to obtain synchronous 
and subsynchronous whirling mot.ons of a horizontal Jeffcott rotor with bearing clearances. The 

Jacobian form for the iterative process which guarantees convergence at all parameter values. The method is showr n 
constitute a robust and accurate numerical scheme for the analysis of two dimensional nonlinear rowr proWe^The 
stability analysis of the steady-state motions is obtained using perturbed equations about the periodic motions. The Floque 
multipliers of the associated Monodromy matrix are determined using a new discrete HB /AFT method. Flip bi ur 
boundaries were obtained which facilitated detection of possible rotor chaotic (irregular) motion as parameters o 
system are changed Quasi-periodic motion is also shown to occur as a result of a secondary Hopf bifurcation due to 
increase of the destabilizing cross-coupling stiffness coefficients in the rotor model. 


Key words: Nonlinear, rotor, clearance, chaos. 


1. Introduction 

Many rotor dynamic systems exhibit nonlinear behavior due to bearing clearances, squeeze film 
dampers, seals and fluid dynamics effects. Nonlinear rotor systems involving bearing clearances 
were studied by several investigators. Bently [1] used a simple horizontal rotor model with a 
bearing clearance to explain the occurrence of subharmonics in his experimental results. Childs [2] 
used a perturbation technique to study the occurrence of subharmonics, assuming small non- 
linearity for the bearing clearance. Saito [3] utilized a harmonic balance method (HBM) along 
with a fast Fourier transform (FFT) procedure, which was originally used by Yamauchi [4], to 
explain some nonlinear characteristics in a Jeffcott rotor on nonlinear supports. Choi and Noah [5] 
also used the HBM with FFT to show the occurrence of super and subharmonics in a rotor in 
presence of a bearing clearance. In [3] and [5], numerical differentiation was used within each 
iterative cycle. This did frequently lead to difficulties in getting consistent convergence in all 
parameter ranges. A numerical approach based on a collocation technique was adopted by 
Nataraj and Nelson [6] and used to obtain periodic whirling motions in nonlinear rotor systems. In 
their approach, the calculation of eigenvalues and eigenvectors is required to obtain steady-state 
rotor whirling motions. This could have the disadvantage of making the numerical process more 
elaborate and lengthy. Nevertheless, the method appears to be versatile and effective. Ehnch [ J 
used numerical integration to show the occurrence of higher subharmonics (up to 9th order) in a 

high speed rotor system with a bearing clearance. 

Simulations revealing aperiodic whirling motion were reported by Childs [8]. Day [9] 
proposed an interpretation involving a ‘nonlinear natural frequency’ to explain the occurrence of 
aperiodic motion obtained using the multiple scales method. 

Nonlinear Dynamics 1: 221-241, 1990 

© 1990 Kluwer Academic Publishers. Printed in the Netherlands. 
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Few analysts have addressed the stability of periodic or subharmonic responses of nonlinear 
rotor systems despite its considerable significance in the development and analysis of modern high 
performance rotor systems. Most of the stability or bifurcation analyses were concentrated on one 
dimensional problems where motion patterns are assumed apriori (Shaw and Holmes [10] and 
Natsiavas [11]). These approaches could be proved unfeasible to extend to two dimensional 
nonlinear rotor problems in which, say, whirling motion involving intermittent contact with a 
bearing clearance would occur. 

This paper addresses the response and stability of a modified Jeffcott rotor system with a 
discontinuous nonlinearity (bearing clearances). The paper consists of two parts. First, a modified 
HBM is developed which combines an exact Jacobian matrix and a Galerkin procedure to 
formulate a robust iterative procedure for determining the periodic solutions. Second, a new 
approach for the stability analysis of the periodic whirling is developed and applied to conduct 
bifurcation analysis of the rotor system and search for possible chaotic responses. 


Equations of Motion 

The equations of motion for a horizontal Jeffcott rotor with bearing clearances (refer to Figure 1) 
can be written as 

mX- + cX' + k,X + Q, Y + x( 1 - y==) " ***** H ‘ “ ^f=! ) 

= meo)~ cos a )t , (la) 




Fig. I. Jeffcott rotor model with bearing clearances (refer to [7]). 
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mY" + cY' + k s Y - Q,X + < Pk h y(l - y==p ) + F^ h x(\ " =2 ) 

2 ■ , (lb) 

= meu) sin a>t - mg , 

where A, is the shaft stiffness, Q s is the cross coupling stiffness, c is the system damping, /x is the 
friction coefficient. 8 is the radial clearance of the bearing. A prime represents a derivative with 

respect to time t and 

. fi, V* 2 + y 2 >g , 

* [o, Va' 2 + y 2 « 8 . 

To study the effect of the parameters on the behavior of the system, the following 
nondimensional groups are introduced: co n = V JChn, K = 4k s k b /(\fk~ s + Xle, y - Yle , 

ft = w /a>„, £ = c/2mw n , y= QJK, a = k b /k s , 8* = 8/e, <}> = g/<o 2 n e, r* = \jx 2 + y\ and vO = coL 
Here v represents the subharmonic ratios. (v=l for harmonic and superharmonic cases, and 
v = n for an nth subharmonic case.) Equation (1) can now be written as 


X + ~ x + £ X + y ^ y + 7X0) - pF{6) = v 2 cos vG , 

ft ft' 4a ft 


2 £v . 

y + ~n y+ n i 


v (1 + v ^) v — y ~ x + F(8 ) + /x T(G) — v 2 sin vd — <t> 2 » 

4a ' ft* “ 


(2a) 

(2b) 


where a dot represents a derivative with respect to the nondimensional time G and is unity if r * 
is greater than 8*, otherwise it is zero. T{8) and F(G) are given by the following expressions. 


y 2 (1+V5)‘ 


xll- 




v 2 (1 + V5) 2 


y 1 - 


yfx 2 + y~ 

8 * 


\f? + 


y 


where 




1, yfx 2 + y 2 > 8* , 

0, VrT ~y-^8*. 


After reaching the steady-state, and assuming periodic whirl, the solution forms of x and y can be 
represented as 


s 


x(G) = a x0 + £ (fl„ cos nG - b xn sin nO ) . 

n = 1 

N 

y(0) = a 0 + 2 («„ cos nG - b yn sin nG ) . 

n = l 

The nonlinear restoring forces of T(0) and F(G) are also expressed as 


(3a) 

(3b) 
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7(0) = c l(t + E (c rn cos n6 - d xn sin nd ) , (4a) 

n = I 
.V 

7(0) = c 0 + E (c yn cos n6 - */„ sin n0) . (4b) 

n = 1 

In equations (3) and (4), N represents the maximum number of harmonic terms considered. 
Inserting equations (3) and (4) into (2), and equating the coefficients of sin(n0) and cos(n0) on 
both sides of the equations, one arrives at the following implicitly nonlinear algebraic equations 


for the constant series terms 

.. . v~ (1 + Va) v , A 

#0) — ^2 4^ a x 0 + 7 a y 0 + C x 0 t LC y(\ ~ ® 

... V 2 (1 + Va) 2 V 2 , . , J. v -n 

g ( 2 ) = 4^ fl v0 - 7 ^2 a x0 + C y0 + nc x0 + 

for the trigonometric series terms 


(5) 

(6) 


2 fvn 


g(4n - 1) = -na xn - -jp b xn + 


v 2 (1 + Va) 2 


4a 


a xn + y 773 a vn + C xn ~ V C yn ~ ¥(*) ^ = 0 

(7) 


, 2 (vn v 2 (1 + Va) 2 

*< 4 ") ' " b - ~ “ft- 4a 


7 ^ b ~~ d 4- urf = 0 

A ^2 jc/i “ y/i 


2£yn v 2 (1 + Va)~ 

IT b ’- + n 5 


g(4n + 1) = -n\„ - -7^ + 77 ^-77 ~ a yn - y ^ a xn + c yn + fic xn = 0 


( 8 ) 

(9) 


g(4n + 2) = nb vn - 


2 fvn 

n 


4a 

° vn + 7? b y " + y 7r- b * n ~ ~ _ 'W 1 ' 2 = 0 • 

“ “ (io) 


In the above equations, 'P(n) is unity if n = v, otherwise 'P(zi) has zero value, and n 
Let the unknown vector P of the displacement coefficients be defined as 

P = [fl x 0’ fl v0’ a x\' b xl ' a yl * b y\ ' • • • i °yN' ^vv] 

and the unknown restoring vector Q of the force coefficients be expressed as 

Q — [CjfO’ Vo * Vl » ^ xl ’ Vl > ^vl * • • • ’ Va/’ .V ] j 


1,2 N. 

(11a) 


(lib) 


where 7 stands for the transpose. The Newton-Raphson method can be used for this two- 
dimensional rotor problem to solve for the unknown vectors P and Q, Alternatively, using 
equations (5 )— ( 10) another iterative scheme such as the Broyden method [12] can be used to 
obtain the steady-state solutions in which calculation of the Jacobian matrix would be avoided. 
Broyden method converges more slowly (usually it requires more iteration steps) but possesses 
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larger radii of convergence for initial guesses. In this study, Newton- Raphson method is used, 
since an explicit form of the Jacobian was made available. 


Newton-Raphson Approach 

Equations (5)-(10) are nonlinear algebraic equations whose solutions yield P. A system of linear 
equations for the correction increments AP of the unknown coefficients can be written as 


[7]AP + G = 0 , 


( 12 ) 


where [7] = [dG/dP] is the associated Jacobian (matrix of first order derivatives) whose elements 

are listed in Appendix A, and G is a (4 N + 2) column vector whose element g(l), g(4n + 2) 

are given by equations (5 )— ( 10) . 

Using an AFT method [13], the nonlinear force vector Q can readily be obtained from the 
unknown vector P. An IDFT is first employed to obtain discrete displacements of x and y from P 
which in turn are used to calculate corresponding discrete values of the nonlinear forces. A DFT 
procedure is then used to calculate the Q vector from these discrete nonlinear forces. As Q is a 
function of P, the Jacobian matrix, [/], has the components of dQ/dP which are expressed as 


1 , M-\ 

dc - = J_ y A 

M A ' 

I M - 1 

- - M 2 


da 


dc 


2-rrlr 2i rnr 

cos . - cos 


M 

2-rrlr 

COS — 7 T- cos 


dd. 


r = 0 
M- 


da. 


dc i. 


j yi - i 

= ^ A, 

iVI r = 0 


COS 


M 

2-rrlr 


sin 


M ’ 

2-nnr 

M 

2-rrnr 


1 2-rrlr . 2-jrnr 

— 2j B r cos ~ sin 


da vl M r=0 


da xl M ' 




ac„ 


^ M - 1 

da vl = M r ? 0 ° r 

dd vn _ J_ *y' _ 

da xl M ~ o ' 

dd - 1 n 

— — = --Tj 2 j D r C os 

da yl M r=0 

«, / = 1, 2, . . . , N, where 

(1 + Vaf 


M 
27 rlr 

cos — — cos 


M ’ 

27 rnr 


M 


lirlr 

cos — cos 
M 


lirlr . 
cos — — sin 
M 

lirlr 


M 

lirnr 
M ’ 

lirnr 
M ’ 

lirnr 


M 


sin 


M 


dc 

1 

M - 3 

V' 



2-rrlr 


2-rrnr 

xn 


Z 

A, 

sin 

M 

cos 

— 

db tl ~ 

~M 

r = 0 



M 

dc 

1 

M- 1 



2-rrlr 


2-rrnr 

xn 


Z 

B 

sin 

M 

cos 


’ db„ " 

= M 

r = 0 

r 


M 

dd rn 

1 

M - 1 
V' 



2-rrlr 


lirnr 

xn 


Z 

A, 

r sin 

M 

sin 

— 

’ db x , ‘ 

" M 

r = 0 

i 



M 

dd Tn 

1 

M- 1 

B r 


2-rrlr 


2Ttnr 

xn — 


Z 

sin 

M 

sin 

— 

db yl ~ 


r = 0 

r 


M 

dc y n 

1 

M- 1 
V 



2-rrlr 


lirnr 



Zr 

C, 

sin 

_ . cos 


db xl " 

= M 

r = 0 


M 


M 

dc yn 

1 

M- ! 



2rrlr 


lirnr 



Z 

D, 

sin 

M 

cos 

— 

dbyl ' 


r ~0 

r 


M 

dd yn 

1 

M- 1 
V' 



2-rrlr 


2-rrnr 



z 

c. 

sin 

M 

sin 

— 

db x , ~ 

~M 

r = 0 

r 


M 

dd vn 

1 

M- 1 

VT 



2-rrlr 


2-rrnr 



z 

D, 

, sin 

M 

sin 

* 

’ db Yl ■ 

" M 

r = 0 

i 



M 


( 13 ) 


[ v- (1 + Va) 2 v 2 (1 + Va) : 2 2U-V2) 2 

^ — w 4 a ” 6 (x +> } y V 
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B, 

C, 

D 


v_ (1 + Va) 
ft 2 4a 

■ - [_Zl (1 + v'a) 

r L n 2 


- 8*(x^ + y 2 ) ( - y2) xy 


8*(x - + y 2 y- 32) xy 


4 a 


,-fr 


i l t _ Zl 11 + ^) : .u ..2M-3/2) 2 


4a 


ft 2 4a 


«*(x‘ + /) 


i 


and M is the total number of discrete data points in the time domain. More details about the 
calculation procedure can be found in Appendix B. 

The procedure of using the Newton-Raphson method to determine a periodic solution can be 
summed up as follows: 

(1) Assume an initial value, P (0) , of the coefficient vector P. 

(2) At a given iteration step, evaluate Q ( *’ from P U) by using the AFT method. 

(3) Calculate [/] and G. 

(4) Solve equation (12) to determine the correction vector AP. 

(5) End iteration if (AP U) — AP ( * _,) ) is within a specified error bound, otherwise set P a + I) = 
P U) + AP </c) and return to step (2). For obtaining possible multiple solutions, different 
initial guesses could be selected at step ( 1 ). 


Stability Analysis 


One of the advantages of the HBM with an AFT procedure is that it readily provides stability 
criteria as well as information concerning bifurcation behavior. In rotor systems, stability and 
bifurcation analysis of a given periodic solution can offer valuable design inputs to avoid sudden 
change of behavior, irregular (chaotic) motion, and dangerous subsynchronous or supersynchron- 
ous vibrations. To investigate the stability behavior of a 27r -periodic solution, eigenvalues of the 
associated monodromy matrix are utilized [14]. 

For the stability analysis, the second order nonlinear ordinary differential equations of the 
present two dimensional problem are perturbed about the determined periodic solution under 
consideration. This leads to the following perturbed equations 


2£v . v 

Ax + — — Ax + — ^ 
ft fl- 


2£v . v 2 

A V + ^ Ay + ^ 


(1 + Vaf_ 

4a 

(lfVi ) 2 

4a 


Ax + y ^3 Ay + A Ax - BAy + fiCAx - fiDAy = 0 , 

1 L 

-> 

Ay - y ~~2 Ax — CAx + DAy + fiAAx - fiBAy = 0 , 


(14a) 

(14b) 


where A, B, C, and D have the same expressions as given previously. Equations (14) are ordinary 
differential equations with periodic coefficients, since A, B, C, and D are 2tt - periodic. Equations 
(14a) and (14b) are cast in first order form, or 

P=["(0)]P, (15) 

where [Ax, Ay, Ax, Ay] r , and [u( 0 )] is the matrix defined as 
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° \ 

1 

0 

~Uv 

n I 


and 

< 7 i( 0 ) 


v 2 (1 + V5)‘ 


fL 


4 a 


+ A + [iC 


/° 

0 


[«(*)] = 


~<l i( 0 ) 

~<h{Q) 


o 

0 


1 

0 


a 


-q*W) 0 


q 2 ( e ) = 7 B ~^ D 

2 

q 3 (0) = -y “2 ~ C + fxA 


1 



(1 + \/a ) 2 
4 a 


+ D- fiB . 


(16) 


Let the monodromy matrix be denoted by [/?], and satisfy the following ordinary differential 
matrix equation 

[*] = MW*]; [*(0)] = (/] . 

where [/] is the identity matrix. Without loss of generality, the initial conditions are assumed as 
the identity matrix. The monodromy matrix can be calculated by integrating equation (17) 
numerically from time 0 to one period. 2tt. The eigenvalues of the monodromy matrix are the 
Floquet multipliers which are used to determine the stability of the 27r-periodic solutions as 
follows, [15], 

1. If all the multipliers are located within the unit circle, the system is stable. 

2. If one of the multipliers leave the unit circle through -1, this indicates period multiplying 
bifurcations. 

3. If one of the multipliers leaves the unit circle through +1, this could indicate bifurcations, 
possibly including a saddle node. 

4. If a pair of complex conjugate multipliers is leaving the unit circle, a Hopf, or a secondary 
Hopf bifurcation could occur. 


Numerical Results and Discussion 

Among the seven nondimensional parameters (fi, £, y. a, 8 *, <f>, /i). the magnitudes of 8* = 30 
and <f> = 30 x stiffness (=(1 + Vaf/Aa) were selected so as to satisfy the condition that the rotor 
center offset equals to the clearance (normal tightening condition [2]). The other five parameters 
were varied. The normal tightening condition not only reduces the number of parameter variation 
effects to be studied, but also fulfills the same whirling motions which were reported experimental- 
ly. This condition is necessary for intermittent rotor/bearing contacts to occur, constituting the 
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main nonlinearity of the system in the y direction. Figures 2, 3, and 4 display the same whirling 
shapes as obtained by Ehrich [7] within the same parameter ranges, as will be discussed below. 

Periodic Response 

The accuracy of the HBM/AFT utilizing a Newton-Raphson algorithm (hereafter the HBM/AFT 
is used to indicate the HBM/AFT with Newton-Raphson for convenience) is compared with 
numerical integration (4th order Runge-Kutta) as shown in Figures 2, 3 and 4. Figure 2 shows a 
period- 1 whirling orbit at ft = 1.1. The figure shows very good accuracy of the HBM/AFT. Figures 
3, 4 show a period-2 (2nd subharmonic) and a period-3 (3rd subharmonic) whirling response at 
Cl = 2.2 and Ti = 3.2, respectively. Again, these figures show the HBM/AFT method to be very 
accurate. Note that small discrepancies in higher subharmonic orbits are due to truncation of 
higher harmonic terms in the assumed steady state solutions. For the results presented herein, up 
to 4 harmonic terms were considered which combined good accuracy with high computational 
efficiency. The other iterative scheme of Broyden also converges to the same orbits as shown in 
Figures 2, 3, and 4 with comparable accuracy. The major difference between these two methods is 
that the HBM/AFT converges much faster than the Broyden but requires more narrow domain of 



X 


Fig. 2. Orbit-1 whirling motion (o =25, £ = 0.02, ft = 1.1, y = 0, /x = 0) — HBM; . . . Runge-Kutta. 
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Fig. 3. Orbit-2 whirling motion (o = 25, ( = 0.02, il = 2.2, y = 0, ft - 0) — HBM; . . . Runge-Kutta. 


initial guesses and more complicated formulations involving Jacobian calculation. However, with 
the HBM /AFT previously calculated results can be used to guess next initial starts for consecutive 
calculation. It was therefore concluded that the HBM /AFT method constitutes more effective 
means of obtaining bifurcation boundaries. 

Bifurcation Behavior 

One of the major advantages of implementing the HBM/AFT method is that it can readily lead to 
a procedure which yields stability and bifurcation boundaries at which qualitative changes in rotor 
whirling occur. 

First, effects of the magnitudes of the stiffness ratio a and critical damping £ were 
investigated. The results show that an increase in a causes period doubling through flip 
bifurcation. Boundaries between stable period-1 whirling motion and stable period-2 orbits are 
shown in Figure 5. In this figure, a stable period-1 orbit exists outside of each curve and period-2 
orbits exist inside of each curve. This figure also reveals that higher £ may eliminate dangerous 
period-2 orbits with the same frequency. This result well agree with previous results [16], 

Figure 6 shows the same a and £ influence on flip bifurcation with ft = 1 . 6-3.0. The figure 
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Fig. 6(a). First flip bifurcation boundaries in a - H plane (Cl- 1.6-3. 0, y - 0, ^ -0). 



Fig. 6(b). First flip bifurcation boundaries in a-H plane (0-2.55-3.0, y-0, Pj period-1 whirling, 

P, = period-2 whirling. 

reveals that there are two types of period-2 orbits possible in the range of ft = 2. 5-3.0, since there 
are two flip boundary branches with fixed Cl. Next, the maximum magnitudes of Floquet 
multipliers are calculated for Cl = 2.7 for different values of a and £ as shown in Figure 7. In this 
figure, there are two types of period-2 orbits (denoted as type A and type B) which are possible 
with £ less than 0.1. These two types of period-2 whirling motions are confirmed by numerical 
integration as shown in Figure 8. Type A response could be considered to be more dangerous 
since it has larger amplitude. Further increase of a leads to another flip bifurcation (2nd flip 
bifurcation) as shown in Figure 9. This figure shows a similar £ effect as that observed in Figure 5. 
In this figure period-4 orbit exists outside of each curve and period-2 orbit is located inside of each 
curve. It is interesting to note that at the range of ft = 1.8-2. 2, higher subharmonics are difficult 
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to obtain unless a has very high value, which approaches an impact condition. From Figure 9 it is 
predicted that by increasing a, further sequences of period doubling occur leading to irregular (or 
chaotic) whirline motion of the nonlinear rotor system studied herein. Figures 10 (a)-(c) show this 
period doubling" process at fl = 1.6. Figure 10 (d) shows chaotic whirling with a high a value. This 
chaotic motion is quite different from aperiodic whirling motion (which is discussed later). The 
occurrence of both types of motion is confirmed by stroboscopic snap plots at every forcing 
period, which is similar to the Poincare maps in one dimensional problems. 

The important characteristics of chaotic motion in the present rotor system are associated 
with its violent vibration which might cause severe rotor-stator interaction. Chaotic motion is also 
characterized by a wide-band, continuous frequency content which might lead to adverse 
conditions of fatigue or excitation of other coupled structures to the rotor. A remedy of this 
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X 


Fig. 10(c). Orbit-4 whirling motion (f) = 1.6, a = 50, £ = 0. 1 , y = 0, fi = 0). 



X 

Fig. 10(d). Chaotic whirling motion (Q = 1.6, a = 100, £ = 0.1, y = 0, ix. - 0). 


situation could be to increase the critical damping or to decrease the shaft-to-support stiffness 
ratio. 

The effect of the friction coefficient, /a, between rotor and stator, is investigated and the 
results are shown in Figure 11 for 0 = 1.5. The figure shows that higher /i tends to stabilize 
whirling near the flip bifurcation boundary. However, it is apparent that y. has little effect on the 
whirling magnitude or stability within stable orbit regions. Figure 12 shows a critical example 
revealing how /i affects the whirling motion near the first flip bifurcation boundary region. The 
figure shows that by increasing /u, the period-2 orbit becomes period-1 orbit but the whirling 
amplitude does not change. Therefore, in critical situations, subharmonic vibration could be 
eliminated by increasing the magnitude of fi. 
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Fig. 13. Hopf bifurcation boundaries (o = 1, fi =0). 


Finally the effect of the cross coupling stiffness, y, is investigated and the results are shown in 
Figure 13. It is seen that the change in y results in a different type of bifurcation. A Hopf 
bifurcation can exist in this case (two complex conjugate multipliers leave the unit circle while the 
other two remain inside of the unit circle). In Figure 13, the period-1 orbit exists below each line 
and a Hopf bifurcation occurs above that line. A Hopf bifurcation produces aperiodic (or 
quasi-periodic) motion as shown in Figure 14. The figure shows that the aperiodic motion has two 
different frequency components (which are incommensurate) and much larger whirling amplitude 
than the period- 1 orbit. 



Fig. 14(a). Aperiodic whirling motion due to Hopf bifurcation (o — 1, M ~ 0) y - 0.39; y — 0.40. 
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1.4-1 

i 



Nondimesional frequency 

Fig. 14(b). Power spectrum of aperiodic whirling motion (a = 1, ft =0. y = 0.40). 


Conclusion 


A robust iterative numerical procedure based on the HBM/AFT method has been presented for 
obtaining the periodic responses of a rotor system on nonlinear supports. Modern bifurcation 
theory is utilized to characterize the dynamic behavior of the system. A bifurcation analysis 
method is developed which provides boundaries of parameter regions at which rotor whirling 
change its shape rapidly, resulting in the occurrence of subharmonic, aperiodic or possible chaotic 
motion. 

The results of this study lead to the following observations concerning the dynamic behavior 
of the nonlinear, modified Jeffcott rotor model considered herein as function of its dimensionless 
parameters: 

1. Increasing the bearing to shaft stiffness ratio, a, increases the degree of nonlinearity which 
makes it possible for a flip bifurcation to occur, possibly producing a sequence of period 
doubling motions. 

2. For the same a and Cl, an increase in £ leads to elimination of the subharmonic motion. 

3. With high values of a, occurrence of chaotic whirling motion is possible. This follows from 

( 1 ). 

4. For the parameters considered herein, the coefficient of friction, /Lt, has little effect on the 
subharmonic response. However, higher /x could eliminate the subharmonics near existing 
flip bifurcation boundaries. 

5. Increasing the cross coupling coefficient, y, could cause a Hopf bifurcation to occur which 
may lead to aperiodic whirling. A more systematic investigation of the quasi-periodic 
response of nonlinear rotor systems is needed. A nonzero value of y is necessary for the 
occurrence of aperiodic solution. This is since in this case a limit cycle can exist in absence of 
imbalance forces. A quasi-periodic response then occurs in presence of an imbalance force 
involving a frequency related to that of the limit cycle and the forcing frequency, or 
rotational speed of the shaft. 
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Appendix \ 

Elements of the Jacobian matrix, [ 7 ] 
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dC xn *>Cyn 


db 


' An, An — 1 


yn yn 

*8* n _ nf dd *n 

^Z~~ nt2 ~^Z 


I - d J± g _ 2 dd x 

4nAn db. m 1 db r 


db 

+ ^ 
+ m 


yn 

dd„ 


da. 


dd 


yn 


db. 


r 

1 

c 

’<r 

<*o 

dd xn 


dd y n 

• / 4 n, 4 n+ 1 

da 

da 


da vn 


yn 

y n 


yn 


dg An 


dd xn 

* s<i >- 

^ An, An 4-2 

Myn 

h - 

db yn 

M a**. 

*^ 4 n + 1 , 4 n - 1 

_ ^&An 4 1 

= — , 

f 3 + 

d Sn . 3 C„ 

h U 

a*™ da xn 


j fyUg+l _ . 

J *n + \An~ db 


dc 


yn 


db. 


+ M 


dc. 


db 


J _ d&4n + l = - n 2 + i + 

*'4/1+ 1 ,4n + 1 « T J, T 


xn 

dc „ 


*^4n 4- 1 ,4n + 2 
An+2,4n-l 


dj>4n + 1 
dg4n + 2 




+ M 


dc: 


dc 

= ■"' 2 + 56 


yn 


yn 


+ M 


da 


y« 


dc. 


yn 


db 


yn 


dd 


yn 


/ _ ^^4n + 2 _ 

* / 4n + 2.4n ^ 3 


dd 


yn 


db. 


~ M 


db. 
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^ 4n + 2.4n * 1 


*g. 


An •> 


8d„ 


dd r 


da VH 


= -nU- — M 

5fl v» da > 


J _ d&4 „ 4- 2 _ 2 _ . _ 

•^n + ’.An-: — • n 


db 


yn 


dd vn dd 

vn in 


at 


vn 


a/?. 


where 


_ V (1 + Va) 
1 1 „ ^ 


1 ft : 4a 


2 fi> / 

'* = 1T’ ' 3=y ^- 


n* 


Appendix B 

Calculation of elements of [7] 
Using 


T(x. >•) = <!> ^3 
Fix, y) = 4- ^ 


(1 + Va) 2 
4 





the incremental form is expressed as 

( dT ST \ 

— Ax + — Ay] = AAx - BAy 

A F(x, y) = Ax + — A y j = - CAx + DAy 


(Bl) 

(B2) 


(B3) 

(B4) 


where 


— £{ 
—Si 


(1W3); _ (1 + Vaf 2 + 2 ,» ,} 
4 4 J 


(1 --^ 


,M(l+J^_ (l+vV 8 . (xi Y 3, V ] 

4 4 ) 


Also, the Ax and Ay are 

dx ^ ( Sx Sx \ 

^ X= da~ 0 A ° z0 + ^ V TfT Afl ™ ” IT - ) = Aa *o + 2 (Aa,„ cos nO - Ab xn sin nd ) , 


n = 1 \da 


db 


n - 1 


(B5) 
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Ay = 



ay 


Aft. 


X -v 

) = Afl l() + S (Aa v „ cos nd - A b yn sin nd) . 

' n = 1 


(B6) 


Similarly, from equation (4), AT and AF can be expressed as 

J V 

AT = A c x0 + 2 (Ac,,, cos nd - A d xn sin nd ) , (B7) 

n - l 

A' 

AF = Ac v0 + 2 (Ac vn cos nd - A d yn sin nd ) , (B8) 

n = 1 

From equations (B3), (B4), (B7) and (B8), and using Galerkin’s method, one can get the 
following expressions for determining the elements of [/]. (The utilization of Galerkin’s method 
rather than DFT and IDFT makes it much easier to obtain the SQ/dP for the present two 
dimensional system.) 


f A [Aa x0 + X (Aa xn cos nd - A b xn sin nd) |{cos d, . . . , sin d} T dd 
Jo l n = 1 J 

- f b{a<j 0 + 2 (Aa vn cos nd - A b yn sin ai0)]{cos 0 , . . . , sin nd} T dd = 
A\ A c x0 + 2 (Ac x „ cos nd - A d xn sin nd) |{cos d, . . . , sin nd) T dd , 

- \ c\ka x0 + X (Afl ln cos nd - A b xn sin n0)]{cos 6, . . . , sin nd} T dO 

JO l n = I J 

+ ( d|aa 0 + 2 (Aa cos «0 - A b vn sin n0) |{cos 0, . . . , sin nd} 7 dd = 
J° l n=i J 

[ /i(ac v0 + 2 (Ac vn cos nd - Arf yn sin n0)]{cos 0, . . . , sin nd} 7 dd , 

J° l rt=rl ' J 


(B9) 


(BIO) 


where the upper limit of integration, 7, is 2tt . Using equations (B9), (BIO), the first derivatives, 
0Q/3P, for the Jacobian matrix are obtained as listed in equations (13) in text. 
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ABSTRACT 

A new HB (harmonic Balance)/ AFT (Alternating Frequency Time) method 
is further developed to obtain synchronous and subsynchronous whirling response 
of nonlinear MDOF rotor systems. Using the HBM, the nonlinear differential 
equations of a rotor system can be transformed to algebraic equations with unknown 
harmonic coefficients. A technique is applied to reduce the algebraic equations to 
only those of the nonlinear coordinates. Stability analysis of the periodic solutions 
is performed via perturbation of the solutions. To further reduce the computational 
time for the stability analysis, the reduced system parameters (mass, damping, and 
stiffness) are calculated in terms of the already known harmonic coefficients. For 
illustration, a simple MDOF rotor system with a piecewuse-linear bearing clearance 
is used to demonstrate the accuracy of the calculated steady-state solutions and 
their bifurcation boundaries. Employing ideas from modern dynamics theory, the 
example MDOF nonlinear rotor system is shown to exhibit subsynchronous, quasi- 
periodic and chaotic whirling motions. 

INTRODUCTION 

There has recently been a tendency to increase the power and efficiency of 
rotating machinery. Smooth running of such machinery is often of great importance 
both for mechanical reliability and for user satisfaction. Consequently various 
rotordynamic effects, which in some cases may be due to existing nonlinearities, 
become increasingly important in the design and operation of such machinery. 

A nonlinearity such as due to bearing clearances or rotor/stator rubs may sig- 
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nificantly alter the vibrational response. In the past, some aspects of the associated 
phenomena have been considered for simple rotor systems [1,2,3], experimentally 
[4] or analytically [5,6,7]. Day [2] and Neilson and Barr [3] showed possible occur- 
rence of quasi-periodic whirling of rotors in presence of bearing clearances. Other 
investigators [1,4, 5, 6, 7] also showed subsynchronous whirling motion occurring with 
nonlinear rotors. Very few studies were performed however for multi- degree of free- 
dom (MDOF) nonlinear rotor systems. 

Nataraj and Nelson [8] adopted a collocation approach and the Guyan reduction 
technique to obtain the steady-state whirling motion of a MDOF rotor system 
with squeeze film dampers. Despite of the complex calculation involved and 
lack of anj' stability analysis, their method is one of few existing approaches for 
determining the steady-state motion of MDOF nonlinear rotor systems. Ehrich 
[9] demonstrated experimentally that “chaotic” motion does occur in a high speed 
MDOF turbomachinery. 

The present stud}' consists of two main parts : i) First, the new har- 

monic balance method (HBM), with alternating DFT(Discrete Fourier Trans- 
form)/IDFT(Inverse Discrete Fourier Transform), is adapted for obtaining the pe- 
riodic steady-state whirling motion of a MDOF rotor with bearing clearances. This 
approach has the advantage of offering robust convergent solution algorithms during 
the iteration process. Stability analysis using the present HBM can yield parame- 
ter ranges for which quasi-periodic whirling motion would occur, ii) Secondly, the 
stability analysis for nonlinear MDOF rotor systems is performed based on per- 
turbation involving the harmonic coefficients of the periodic solutions. Through 
application of modern dynamic theory, all possible parameter ranges for the sud- 
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den change of whirling motion (i.e. bifurcation) can be obtained. This stability 
approach cam be extended for any MDOF rotor systems with piecewise-smooth 
or polynomial-type nonlinearities. Another advantage of this approach is that it 
enables predicting parameter ranges in which chaotic rotor whirling would occur. 

SYSTEM REDUCTION 

A typical multi-disk rotor system with nonlinear bearings is considered. For 
a finite element (FE) formulation of the nonlinear rotor system, the rotor shaft 
segments are modeled as Euler beam elements taking shaft rotation effects into 
account [10]. An assumed axisymmetric geometry of the rotor shaft elements leads 
to the same mass and stiffness matrices in the X-Y and X-Z planes. If internal 
damping is neglected, the system equations of motion can be expressed in terms of 
the assembled mass, damping and stiffness matrices [M], [C] and [K] as 

[M]q+ [C]q+ [AT]q = fu + fn (!) 

where q = [t /,z] T and y, z denote the 2Lxl assembled state vectors in the X-Y 
and X-Z planes, where L is the total number of nodes. [C] has nonzero off-diagonal 
elements including damping opposing the motion of the disks and gyroscopic terms 
which couple the motion in the two planes. [Kj has also nonzero off-diagonal 
elements due to rubbing and other coupling forces between rotor and stator. The 
vector f u represents the vector of the disks’ mass imbalances and side forces on 
the disks, while the vector f n denotes the force vector at the nonlinear bearing. In 
particular, the nonlinear restoring forces in the y, z-directions on the bearing with 
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clearances, at the node j, are expressed as 


fy, = hyj{l ~ 


\[v? + 


0 


(2 -a) 


fz, = - 


\/V + 


(2-6) 


where kb is the bearing stiffness and 8 is the clearance between bearing outer race 
and stator. The bearing forces f V} or f z . will vanish if 8 is larger than the radial 
displacement, otherwise the bearing forces will be as given by equation (2). 


In equation (1), the matrix size becomes 4Lx4L, where L is the total number 
of nodes. This may require a very large core size and much computation time 
to calculate the dynamics of the nonlinear rotor-bearing system. In linear rotor 
dynamics, system matrix reduction techniques, such as the component mode 
method [11] or the Guyan reduction technique, could give a reduced form of the 
dyna mi c equation of motion. As an another approach, the impedance method [12] 
can reduce the number of equations of motion to obtain forced responses at specific 
locations of the system from those of “master”, (kept) degrees of freedom. This 
reduction is exact and does not involve any approximation. In a similar fashion, 
large nonlinear rotor dynamical equations can be reduced to obtain steady-state 
response at bearing location using the HBM. This is achieved by using a version of 
an impedance formulation in which the system is reduced to its displacements at the 
bearing clearances, (see [12] and [20]). In the present study, the impedance method 
is applied to each of the harmonic components of the assumed periodic solution. 


Extending the procedure developed by Kim and Noah [14], a periodic solution 
for the motion of the rotor is represented by a finite Fourier series expansion. For 
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the displacement of the i th node (i=l,2,...,L), one can write 


N 


E . nut ■ . nut ^ 

\ a ny COS ^nv _ / 

n=l 


V 


E h , x nUjt ,t • ™ t \ 

i a nz COS — “ Kz Sln — ) 

n= 1 


(3-o) 


(3-6) 


where u is the subharmonic ratio, which is unity for harmonic and superharmonic 
cases, or an integer other than unity for subharmonic cases, and u> is the shaft 
rotational frequency. Since the motion is periodic, the nonlinear restoring force can 
be written as 

(4 - a) 


, J j nU} t ,j . 

hi = c i y + cos ~ ~ d *y sm 


n= 1 
N 


t j i i nu)t . nu)t 

hi = C 0z + cos — - d 3 nz sm — ) 


(4-6) 


71=1 


where j denotes the j th bearing node (j=l,2,...,m). The advantage of introducing 
equation (4) is to avoid the difficulties which would arise if standard harmonic 
balance procedure is directly applied to equation (2). The denominator term of 
equation (2) would not then be easily expressed in harmonic terms resulting from 
using equations (3-a) and (3-b). 


Substituting equations (3) and (4) into equation (1) and equating the coeffi- 
cients of similar harmonic terms lead to the following algebraic relationships : 


constant terms 

[K] A 0 = C° (5) 

where A 0 and C° represent the constant Fourier coefficients of equations (3) 
and (4), respectively, and both are of size 4Lxl. 
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cosine terms 

_(™ 2 )[M]A n -( — )[C]B n + [A r ]A n = C n + F uc (n = 1,2, N) (6) 

where A" and B n represent the n th cosine and sine coefficient vectors, 
respectively, both of 4x1 dimension, C" denotes the n th cosine coefficient 
vector of f n , while F uc represents the cosine coefficient of the imbalance and 

side force vector F u - 


sine terms 

— )[C]A“ - [Jf]A“ = D"+F„. (n=l,2 (7) 

v u n V 

where vectors A.B are the same in equation (6), D n denotes the n ih sine 
coefficient vector of f n , and F u . is the sine term of the imbalance and side force 

vector F u . 


At this stage, one can reduce equations (5), (6) and (7) using the impedance 
reduction technique, only retaining the coordinates at the bearing nodes. For the 
constant terms, equation (5) can be partitioned as 


( [Kkk] 

V [Krk] 



( 8 ) 


where the subscripts “r” and “k” denote reduced out and kept coordinates, respec- 
tively. By applying the elimination procedure, the constant terms can be repre- 
sented in the reduced form 

[K' kh \Al = Cj (S) 


where [.K* fc ] is the reduced stiffness matrix of the rotor and A®, C£ are both 2mxl 
vectors, where m is the total number of bearings. The n tk cosine and sine harmonic 
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terms can be reduced in a similar fashion to yield 

[Sid nA£ = CE + U" (10) 

where [5i t ] n is the reduced matrix involving system mass, damping and stiffness 
matrices, A£ is the n th Fourier coefficient vector of equation (3) at the bearing, C k 
is the Fourier coefficient vector of equation (4) and U is the reduced vector 
involving imbalance and side forces. Here A£, C£, U n are all 4x1 vectors for 

each n. 

Combining equations (9) and (10), the following assembled matrix is obtained 

[T]Ak = Ck + V (11) 

where [T] is a 2m(2N4-l)x2m(2N+l) matrix, A k , and C k are 2m(2N+l) vectors 
which represent the trigonometric Fourier coefficients of equations (3) and (4) at 
the bearing node, respectively, and V is the total sum of the reduced imbalance 
force vector. It is noted here that the total dimension of the system equations is 
reduced from 4L to 4mN-b2m. For a multi-degree of freedom rotor, L is very large in 
comparison with m (number of bearing nodes) and N (number of retained harmonic 
terms), so the system is significantly reduced for obtaining the steady-state response 
of the nonlinear rotor system. Another advantage of the reduction technique using 
the HBM is that stability analysis of the periodic response of multi-degree of freedom 
rotor systems can be made highly efficient by considering only truncated harmonic 
terms. This is discussed in the section further below on stability analysis. 

In equation (11), the only unknown vector is Ak* The vector Ck can be 
calculated using DFT and IDFT relations since Ck is a nonlinear function of Ak* 
This procedure is described next. 


53 



10 


Calculation of Ck from Ak 

To determine C k from A k , it is necessary to first calculate the discrete 
y and z values of yA = [y(A<), y{2At ), ..., y((M + l)At)] T and za = 

[z(A<), z(2At ), ..., z((M + 1 ) At)]" 3 ", where At is any discrete time steps, M+l 
is the total number of discrete points, and superscript T denotes the transpose. 
Using IDFT, the discrete displacements, yA and Za, can be obtained as 





yA = [Va 

l A yk 


(12 - a) 




ZA = [Va 

]A*k 


(12 - b ) 

where 








[Vi] = 

( l 

1 

i 

cos At/u 

0 

— sin At/u 

... 

1 

cos At/u 

1 

cn 

■ 

... »o 

r-U 

(13) 


vi 

cos MAtjv 

— sin M At/u 

cos M At/u 

— sin MAt/u / 



and in which A y k, Auk denote the cosine and sine parts of Ak, respectively. The 
discrete nonlinear restoring force f y A, ft A can be obtained using equations (2) and 
(12) as 

fyA = fy(yA,ZA) (14 - a) 

f«A = f*(yA,ZA) (14-6) 

where 

fyi= [/,(A<),/ > (2Af),...,/ s ((M + l)A!)l T (15 - a) 

f.A = + l)A,)f (15 - b) 

As Ck is the Fourier coefficients of f y A and DFT offers the following expression 

c„ = MTT |Zi|f(ViAk) (16) 
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where 


[Z*] = 

/ 0.5 
1 
1 

0.5 

cos A t j v 
— sin A tjv 

0.5 

cos2Af/i' 
— sin2Af/i' 

0.5 \ 

cos (M -t- 1 ) At / 
-sin(M + l)At/u 


1 

— sin N Atjv 

- sin 21V At/ 1; 

— sin (M + l)NAt/u / 


and f is the sum of f y and f e . 


From equation (11), one can solve the nonlinear algebraic equation by a 
successive iteration procedure. Newton- Raphson scheme is one of the favorable 
techniques for obtaining steady-state solutions since it has rank two convergence. 
The disadvantage of this method is that it requires calculation of an explicit form 
of the Jacobian which is not simply obtainable in multi-degree of freedom systems. 
However, since the problem only involves nonlinear bearing coordinates, the explicit 
form of the Jacobian can be obtained using the DFT and IDFT procedure. One 
of the techniques to avoid calculation of the Jacobian is to utilize the forward 
differentiation [6] which replaces the use of the Jacobian [14]. From the authors 
experience, forward differentiation poses the difficulty of requiring control of the 
differentiation length for each parameter and obtaining the probable convergence 
values near the resonance responses. Moreover, forward differentiation convergence 
rate is much slower than when directly using the explicit Jacobian formulation. 
Therefore, in this study an explicit Jacobian formulation is utilized to enhance the 
computational efficiency and to guarantee convergence for all parameter ranges. 
Another iterative scheme, where calculation of a Jacobian is not required, is to 
use the quasi-Newton method [15]. Although this method avoids calculation of the 
Jacobian, its convergence is much slower when applied to nonlinear multi-degree of 
freedom rotor systems. 
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Newt on- Raphe on Method 

To apply the Newton-Raphson method to determine the unknown vector A k , 
equation (11) can be put in the form 

G = [r]A k -C k -V = 0 (18) 


In using the Newton-Raphson algorithm, the following derivative is needed (using 
equation (16)) 

,c>f(V A A k ) 


dC k 


dA k M + 1^' dA k 
It follows from equations (16) and (18) that the Jacobian column vector of G is 
given explicitly by the relation 


(19) 


[J] = [T] + 


dC k 

dA k 


( 20 ) 


since [J] = and A k is the only unknown vector. The Newton-Raphson 

algorithm for the unknown vector A k with an initial guess A^ ^ can be described 
as 

[7]AA k (p) 4- G (p) = 0 (21) 

where the superscript, p, denotes the p th iteration number. The above algorithm 
terminates after r iterations, so that 


IG 


( r ) I 


< € 


where t is small number, and A^ ^ is a fincil solution. 


STABILITY ANALYSIS 


( 22 ) 
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From the previous section, the steady-state, harmonic or subharmonic whirling 
responses, are obtained utilizing the present HBM formulation. As the steady- 
state whirling responses are calculated using Newton-Raphson technique, a stabilitj 
analysis is necessary to check whether the obtained responses are stable. Multiple 
solutions could exist for a given set of parameters. Some of these solutions could 
become unstable and bifurcate to other forms of solutions. 

If the nonlinear rotor system is of small number of degrees of freedom, it is 
straightforward to perform a stability analysis through perturbation of the periodic 
solution obtained [13]. However, if the nonlinear rotor system possesses large 
number of degrees of freedom, excessive computational time would be required 
to check the stability or to obtain bifurcation information. This is since for the 
stability analysis, the procedure would involve integration of a large matrix in 
order to calculate the monodromy matrix for the perturbed equations with periodic 
coefficients. In the present paper, a more efficient stability analysis method for 
MDOF nonlinear rotor systems is presented. 

From equation (10), the reduced system matrix with the n th harmonic terms 
is rewritten as 

IS1J.AJ = c; + U” (23) 

where [S* fc ] n is the reduced system matrix, A£ is the vector which represents the 
Fourier coefficients of displacement at the nonlinear bearing node, C£ is the vector 
of Fourier coefficients of the nonlinear restoring forces and U n is the reduced vector 
involving the imbalance and side forces. The unknown vectors A£ and C k in 
equation (23), which were already obtained by the HBM method, have the following 
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elements 


A n 

lr — 


cc = 


m m 


L j=l j=l J=1 > =1 

m -i T 


m m 


(24) 

(25) 


E <» • X rf "v ’ S c "- ’ S 

*■>=1 j=i j=i i =1 
where T denotes a transpose, j represents all the j th bearing nodes, and m is 

the total number of bearings. As equation (23) represents the reduced system 


matrix involving only the nodes at the nonlinear bearings, all the coefficients of 
A£ are in general coupled, i.e. the off-diagonal terms of [^kln matrix have non- 
zero elements. At j 1 * 1 nonlinear bearing node, the damping and restoring force are 
affected by all the damping and stiffness coefficients, Cij, Kij, (i,j = 1,2, 
respectively. Therefore, coupled terms in damping and stiffness coefficients can lead 
to the following equations 

m m 

MyjVj + Cyjiifi + kyjiyi — fyj + Uy 

i=l t=l 

m m 

M z jZj + ^ C ZJX z, + ^2 K z jiZi = f.j + Ul (27) 

i=l i=l 

where M, C and K represent reduced mass, damping and stiffness coefficients 


matrix, respectively, fyj , f z j denote nonlinear restoring forces in the y, 2 -direction 
and subscript j represents the j ih nonlinear bearing node. It is noted here that the 
Fourier coefficients of f y and f y are represented by CJJ in equation (23). The next 
step is to calculate each reduced mass, damping and stiffness coefficients from the 
reduced system matrix 


From the previous section, the displacement of the j tfl bearing node can be 


represented by only the n th harmonic component as 

ivt u d 

y J = a hy cosn b ny sinn — 

y u V 


(28 - a) 
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Ujt 


Ijjt 


z 3 = al,cosn b 3 sin n — 

v v 


(28 - a) 

where j = l,2,...,m. Inserting equations (28) into equation (26) and (27), the 


following relations can be arrived at : 

7 /n “' 2 M yl +K yl -C yl *? 


C n u) 

yi~ 


0 

0 


V 


o 

0 


0 

0 


' nw ^2 

C nuj 
zm 


-(?f) 2 M zm +K 


zm 


<y\ 


C ny\ 

b 1 


dl. 

ny 


ny 

„ m 



a nz 


C n: 

\b tJ 


\d?J 


+ u* 


0 

0 

C nu? 
zm v 


(29) 


As equation (23) and (29) are the same, all reduced mass, damping and stiffness 
coefficients of equations (26) and (27) can be calculated from the matrix In 

equation (23). 


Equations (26) and (27) can now be rewritten as 


*-yj 


m ^ m jr 


t=i 


v: 


i=i 


M. 


y> 


/) 


yj 


u v 




M, 


yj 


M, 


= 0 


yj 






+ K zj _ fz± U: 


M S j ^ M.j 

1=1 3 t=i 3 


M; 




M, 


= 0 


(30) 

(31) 


Equations (30) and (31) are the reduced system equations to be utilized for the 
stability analysis of the steady-state whirling motion. It is noted here that the 
system equations are reduced from 2L to 2m where L is the total number of nodes, 
including the number of bearing nodes, m. Therefore, for MDOF nonlinear rotor 
systems, the reduction of system equations has the effect of rendering the stability 
analysis significantly more efficient. 
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Perturbation of equations (30) and (31) can be obtained using Taylor expansion 


as 


dh y j A « dhyj A • L Y^ dhyj 

a^ = ^a w + E^-a w + E^ 


A* + + ^A* = 0 (32) 


3f, 


dyj 


dzj 


^ ^ + E ft* - E ft** + ft* + ft Ai ' * 0 (33) 

J i= 1 *— 1 


where j=l,2,...,m. 


Equations (32) and (33) are transformed into a system of first order equations 

with 

u = [Ay,, Az jt A yj, Ai ; ] T (j = l,2,...,m) (34) 

The resulting variational equations can be written in the following form : 

u = [W(t)J u (35) 


where u e R 4m and W(t) is a 4mx4m matrix. The stability problem of the 
prescribed motion can be formulated as the local stability analysis for u = 0. The 
fundamental matrix \Z{t)\ for the ordinary differential equation (35) with periodic 
coefficients is related to [Z(t + T)], which is also a fundamental matrix by [18] 

[Z(t + T)] = [Z{t)][S] (36) 


where [5] is referred to as the monodromy matrix for the system. The monodromy 
matrix [5] can be produced by evaluating the matrix [Z] at the end of one period 
for the system. The Floquet multipliers, or the eigenvalues, p, of S determine the 
stability of the system. When all the multipliers are of absolute value less than 
unity except for some with absolute value of unity, then the system is located at the 
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stability boundary. Three different types of loss of stability of a periodic solution 
can occur [19] for the system considered when one or more fx are given by 

i) fi=l, saddle-node or transcritical 

ii) p=-l, period multiplying 

iii) fi = ct + z/3 1 = 1, secondary Hopf bifurcation 

The bifurcation can be super- or subcritical. 

AN EXAMPLE MDOF ROTOR SYSTEM 

To demonstrate the application and computational efficiency of the new 
HBM/AFT, the method is applied to a simple MDOF model shown in Figure 1. 
The equations of motion of the rotor are formulated using the finite element method 
(FEM), employing Euler beam elements. The rotor is supported on a linear bearing 
at the left and by a nonlinear bearing with a gap at the middle. The rotor is sub- 
jected to imbalance forces and a constant direction side force. The detailed rotor 
configuration is shown in Table 1. 

First, the whirling orbit at the nonlinear bearing is obtained for a rotor with 
2000 rad/sec rotational speed, as shown in Figure 2. The solid line in this figure 
represents the HBM solution which is seen to be accurate. The minor discrepancies 
between the numerical integration and the HBM solutions are due to : i) truncation 
errors in assuming a finite number of harmonics for the steady-state solution and for 
the nonlinear restoring force. This error can be reduced by retaining larger number 
of harmonic terms, ii) errors introduced by the Guyan reduction. 

A more complicated subsynchronous whirling motion at a rotor spinning speed 
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of 1600 rad/sec is shown in Figure 3. The solid line represents the response obtained 
by the HBM, while the dotted line stands for the solution of direct numerical 
integration. The results show good agreement between the HBM and the numerical 
integration methods. In most of the simulations, only four harmonic terms were 
considered in obtaining the whirling response. This results in drastic reduction in 
the computational time in comparison with direct numerical integration. 

• Another characteristic of the behavior of the nonlinear rotor is the sudden 
change of its whirling shape with small changes of certain parameters. This is due 
to bifurcation whose study is very important since it can lead to a sudden change 
from synchronous to subsynchronous motion (including subharmonic or quasi- 
periodic motion) or vice versa. To better understand the nonlinear MDOF rotor 
characteristics with bearing clearances, stability charts, or bifurcation boundary 
plots, are utilized. Figure 4- a shows the effects of the gap (clearance) and mass 
eccentricity on Hopf bifurcation boundaries. The regions marked by “A" indicates 
quasi-periodic whirling motion. These regions were obtained using present stability 
analysis method which indicated a Hopf bifurcation at their boundaries. Region B 
indicates a stable harmonic whirling motion, i.e. all the roots of the monodromy 
matrix are located within the unit disk. To confirm these stability boundaries, 
numerical integration was performed for a selected set of parameter values, and the 
results are displayed in Figure 4-b and 4-c. These results show that the stability 
analysis developed here is quite accurate. The figure shows that for the range of 
parameters utilized, a complicated quasi-periodic whirling motion can be eliminated 
by decreasing the bearing clearance, while changing the bearing stiffness does not 
affect the motion. 
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Figure 5-a shows effects of the rotational speed and side force on the flip (period 
multiplying) bifurcation boundaries. Region A of this figure represents the primary 
stable whirling motion (i.e. the stable whirling when the rotor does not contact 
the gap), region B is the subsynchronous whirling of order 1/4 which can cause 
violent whirling motion, region C is the subsynchronous whirling with order 1/3, 
and region D shows the secondary stable synchronous whirling motion (i.e. the 
stable whirling when the rotor stays in contact with the supports through the gap). 
Numerical integrations for the orbits in each regions are also performed to confirm 
the bifurcation boundaries as shown in Figures 5-b-e. Since the MDOF rotor system 
displays a 550 rad/sec critical speed with the rotational speeds around 3 or 4 times 
this critical speed, the subsynchronous whirling motions might, therfore, occur. 
Another interesting phenomenon emerging from the results in regions B and C is 
that the motions can bifurcate further by increasing the side force (i.e. by period 
doubling or tripling according to the side force variation). As the HBM can handle 
only finite number of low order subharmonic terms (so that chaotic whirling motion, 
which possesses a large number of low order subharmonics, can not be obtained 
using the HB/AFT method), numerical integration is used to study further flip 
bifurcation in the region C, as shown in Figures 6-a-e. Since region B has similar 
flip bifurcation pattern, the results of this region is not included here. Shown in 
Figure 6-a is the primary stable whirling motion. By increasing the side force, 
the shape of whirling orbits continue to change until subsynchronous motion with 
order 1/3 is achieved. This is shown in Figure 6-b. Further increase of the side force 
results in bifurcation to subsynchronous motions of order 1/9 and 1/27, as shown in 
Figures 6-c and d, respectively. With a small increase of the side force beyond these 
subsynchronous regions, whirling can become almost non-periodic ( chaotic ) as 
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shown in Figure 6-e. Figure 6-f, which shows the response spectra of that in Figure 
6-e, clearly indicates chaotic motion. The reversed bifurcation (i.e. inverse process 
of above bifurcation from primary stable whirling to chaos) begins to occur suddenly 
(“crisis”) by changing from chaotic motion to secondary stable whirling when the 
side force is increased, as shown in Figure 7. 

Figure 8 shows the effects of the mass eccentricity and side force on the flip 
bifurcation boundaries. The region marked by A has a primary stable harmonic 
whirling, while region B has a subsynchronous whirling shape. The figure shows 
that the primary whirling shape does not change its motion very rapidly with 
the eccentricity or with the amount of imbalance. However, the secondary stable 
whirling becomes unstable for the same side force. All bifurcation regions were 
confirmed using numerical integration. The results of the comparisons are not 
included here since they display similar whirling shapes as those of Figure 5. Figure 
9 shows the effects of clearance and the side force on the flip bifurcation boundaries. 
The regions A,B, and C of this figure have the same whirling motions as those in 
Figure 7. The figure shows that larger side force is needed to cause primary or 
secondary stable harmonic whirling when the gap size is increased. 


CONCLUDING REMARKS 

A further developed HBM, using DFT /IDFT, is employed to obtain the steady- 
state periodic response for MDOF rotor systems with bearing clearances (piecewise- 
linear nonlinearity). The Guyan reduction technique is utilized to reduce the system 
to the nonlinear coordinates. The stability analysis is performed via perturbation 
of the obtained periodic solutions. The reduced and approximated system param- 
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eters (mass, damping, and stiffness) are calculated from the determined harmonic 
coefficients. 

A simple MDOF rotor system with a bearing clearance is used for illustration of 
the method. The results obtained show that i) the HBM method, as developed here, 
is robust and efficient, ii) the method leads to accurate bifurcation boundaries for 
nonlinear MDOF rotor systems. Furthermore, the method can in general be applied 
to MDOF rotor systems with piecewise-smooth or polynomial type nonlinearities 
at the bearing supports. 

If the system possesses quasi-periodic response for a given set of parameters, 
the present HB/ AFT approach cannot be used to obtain the corresponding whirling 
motion. However, a modified DFT, which is now under study, is believed to be able 
to produce quasi-periodic whirling motions. 
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Figure 1. MDOF rotor model with a piecewise-linear 
bearing clearance 
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Figure 2. Comparison between HBM and numerical integration 
(speed=2000 rad/sec. gap=3 nun, side force=1400 N. 
eccentricity =1 mm, damping=10 N-sec/m, bearing stiff. = lE+06 N/m) 
HBM ; ... Numerical integration 
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Figure 3. Comparison between HBM and numerical integration 
(speed=1600 rad/sec, gap=3 mm, side force=1400 N, 
centricitv=23 mm, damping=10 N-sec/m, bearing stiff. -1E+06 N/m) 
HRM • ... Numerical integration 




ECCENTRICITY (cm) 


Figure 4. Effect of gap and eccentricity on Hopf bifurcation boundaries 
(speed = 2300 rad/sec, side force=2800 N. dampmg=10 N-m/sec) 

(b) - Region A (gap=3 mm, ecc.=20 mm) ; 

(c) - Region B (gap=3 mm, ecc.=40 mm) 


(a) 
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SIDE FORCE (N) 


Figure* 5. Effect of side force and speed on 
( g ap=3 mm. damping=10 N-m sec. 


flip bifurcation boundaries 

stiff. = lE+06 N/m) 


lb) • Region A (speed=2100 rad/sec. side !.=400 N) ; 
(c) • Region B (speed=2100 rad, sec. side f =500 N) ; 


(d) - Region C (speed=1600 rad/sec, side (.=500 N) ; 
(e) - Region D (speed=2100 rad/sec, side (.=600 N) 
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Figure 6. Flip bifurcation process with speed of 1600 rad/sec 
(a) - primary stable whirling motion (side f.=450 N) ; 

(b) - subsvnchronous whirling motion with order 1/3 (side f.-500 N) 

(c) - subsvnchronous whirling motion with order 1/9 (side f.=52i N) 

(d) - subsynchronous whirling motion with order 1/2 < (side f. 535 N) 

(e) - chaotic whirling motion (side f.=550 N) ; 

(f ) . power spectrum of chaotic motion (side f.=550 N) 


(a) 
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Figure 7. Bifurcation plot due to side force 

+ 1F4-D6 N/m, speed=1600 rad/sec) 

(gap=0-03 cm, damping 10 N-sec,m. s«ft.=lE+06 N/m, 

Pj primary synchronous whirling motion 
Pj subsynchronous whirling motion with order 1/3 
p, subsynchronous whirhng motion with order 1/9 
Pl ’ secondary synchronous whirling motion 
_ QBMf stable solution) , .... HBM(nnstable solution) , o numencal integration 
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Fig „re a. Effect ot side iorce and eccentncity on * «— ^ 
speed=2100 rad/sec, gap=0.03 cm. dampmg=10 N-eec/m, st.ff. 1 + 

\ ,„ m er, whirling , B Asynchronous whirling ; C seconder, whirhng 



ECCENTRICITY (cm) 
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SIDE FORCE (N) 


Figure 9. Effect of side force end gap on flip ^cation boundaries 

, - in N-sec/m, stiff.= lE+06 N/m) 

lspeed=2100 rad/sec, eccentriaty=0.001 cm, ampi 

A ^ whirling ; B subsynchronous whirling ; C secendae, wfurbng 
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APPENDIX C 

"A Convolution Approach for the Transient Analysis of Locally Nonlinear Rotor Systems," 
ASME Journal of Applied Mechanics, Vol. 57, pp. 731-737, 1990. 
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A Convolution Approach for the 
Transient Analysis of Locally 
Nonlinear Rotor Systems 

A computationally efficient convolution method, based on discretized impulse re- 
sponse and transition matrix integral formulations, is developed for the transient 
analysis of complex linear structures interacting through strong local nonlineanties. 
In the formulation , the coupling forces due to the nonlineanties are treated as external 
forces acting on the coupled subsystems. Iteration is utilized to determine their 
magnitudes at each time increment. The method is applied to a generic rotor-housing 
model representing a turbopump of a space shuttle main engine ( SSME ). In that 
model, the local nonlinearity is due to clearances between the rotor bearing outer 
races and the carrier attached to the housing. As compared to the fourth-order 
Runge-Kutta numerical integration methods, the convolution approach proved more 
efficient and robust for the same accuracy requirement. This is due to the closed - 
form formulation of the convolution approach which allows for the use of relatively 
larger time increments and for a reduction in the roundoff errors. 


Introduction 

The use of a direct numerical integration for determining 
the transient response of coupled nonlinear rotor-flexible hous- 
ing systems may require excessive computational time and in- 
volve unacceptable computational roundoff errors. To remedy 
these problems, different procedures have been developed by 
analysts to determine the transient response of large order rotor 
systems. The procedures can be recognized as falling under 
one of two basic approaches: those using physical or modal 

coordinates of the complete system and those using the co- 
ordinates of the individuals components of the system. The 
methods also differ in the numerical integration methods se- 
lected for the analysis. 

Rouch and Kao (1980) employed Guyan (static) reduction 
method to arrive at a reduced size model in terms of the re- 
maining physical coordinates. Accuracy of the results could 
be expected to be acceptable, since the rotor is basically a train 
of mass-stiffness subsystems. Nordmann (1975) attempted to 
minimize the inaccuracy of static condensation by applying 
the static reduction technique to an arbitrarily substructured 
rotor system and then assembling the reduced substructures to 
form a reduced system. The procedure is very laborious and 
no guarantees of accuracy are apparent. 
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Childs (1978) utilized free-interface modes of the various 
system components to represent the assembled turbopumps of 
a space shuttle main engine (SSME). The method, using fourth- 
order Runge-Kutta integration, does not allow for size reduc- 
tion of the housing model while maintaining sufficient modal 
coordinates for accurate representation of the housing. Nelson 
et al. (1982), on the other hand, uses fixed-interface complex 
component modes to assemble a reduced size model. For sys- 
tems with a large number of coupling points among the com- 
ponents, the approach suffers from an introduction of higher 
frequencies resulting from an excessive number of constraints 
imposed at the coupling (or boundary) points. In a transient 
analysis, this will necessitate employing much smaller time 
increments and consequently will lead to excessive computa- 
tional time. 

Only a few analysts have presented techniques for the general 
transient analysis of large nonlinear rotor systems. Adams 
(1980) used a normal mode representation for the rotor in 
terms of its undamped, free symmetric modes and treated 
gyroscopic and nonlinear terms as pseudo-external loads. The 
method presented by Childs (1978) makes use of a similar 
procedure to couple the rotor to its flexible housing. Nelson 
et al. (1982) developed a general computer code for the tran- 
sient analysis of large rotor systems. The user may utilize time- 
step integration in the constrained-rotor (fixed-interface) modal 
space. Again, ail connection points, including those at the 
nonlinearities, must he constrained, leading to the same short- 
comings described previously. 

Spanos et al. (1988) considered linear systems analyzed as 
decoupled subsystems. The accelerations at the interface are 
predicted at the beginning of each time step. The accelerations 
of the interface and interior nodes are then calculated and 
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corrected through iterations. The decoupling makes the mod- 
eling of the assembled system more efficient. Hagedorn (1988) 
proposed a method for the analysis of a large linear subsystem 
coupled to one or more nonlinear subsystems. The dynamics 
of the system is represented by an integral equation of the 
convolution type containing the transfer matrix of the linear 
subsystem, and by the functional relation describing the non- 
linearities. It appears that since no iterative technique was used, 
the method would require a very small integration step. 

A method which proves to be highly efficient in determining 
the transient response of linear systems under specified general 
excitations is to utilize the convolution integral, or the tran- 
sition matrix (Meirovitch, 1980) of the system. Von Pragenau 
(1981) utilized the transition matrix, stating that it offers the 
simplicity of the Euler method without requiring small steps. 
Von Pragenau maintained that for systems with constant coef- 
ficients, the stability and accuracy of the method are acquired 
through the closed-form solution of the transition matrix. 

In the present study, the convolution methods (impulse re- 
sponse and transition integrals) arc shown to be very effective 
when extended for application to linear systems with local 
nonlinearities. The forces at the nonlinear locations are treated 
as the external forces on the systems or subsystems and iteration 
is used at each time increment to determine the magnitude of 
the forces for subsequent increments. 

The technique presented in this paper can also be applied 
in terms of the modes of subsystems calculated separately. 
However, system equations of motion will be formed as in a 
component mode synthesis approach. The response of each 
substructure is solved separately by treating the interface forces 
of the system as external forces to the subsystems. The con- 
volution approach can also be applied to the generalized sec- 
ond-order system equations developed from component mode 
synthesis. In contrast to the convolution approach, if local 
nonlinearities occur at the interface, standard free interface, 
component mode synthesis methods fail. 

The convolution approach is applied to a general rotor- 
housing system with rotor imbalance during startup or shut- 
down. In the present work, eigencoordinates are used to rep- 
resent both housing and rotor. The local noniineanties were 
taken as deadband clearances at the rolling element bearings 
which support the rotor in its housing. The integral formulation 
of the rotor motion is represented by its transition matrix and 


that of the housing by a convolution integral (based on the 
housing’s impulse response). 

The convoluted impulse response can only be applied to a 
system of uncoupled equations while the transition matrix for- 
mulation, in addition, can be applied to coupled equations. 
The transition matrix can therefore be applied to coupled dy- 
namical systems represented by their physical coordinates or, 
in case of rotors, coupled by the gyroscopic terms in otherwise 
decoupled modal representation. Kubomura (1985) used a con- 
volution integration method to achieve dynamic condensation 
of a substructure to its coupling points to other structures. 
Convolution was also used by Tongue and Dowell (1983), and 
Clough and Wilson (1979) to reduce system coordinates to that 
at the nonlinearities. 

Modeling of a Rotor-Housing System 

A representative complex rotor system with a flexible hous- 
ing is shown in Fig. 1. The particular model shown in which 
the present method can be readily applied, represents the high 
pressure oxygen turbopump (HPOTP) of the SSME. The in- 
teraction forces between rotor and housing include various 
seal, impeller, turbine tip clearance, bearing clearance, and 
fluid side forces. 

(a) Rotor. The equations of small transverse motion, 
(Noah, 1986), (S) of the rotor under transient external and 
imbalance forces may be written as 

l+m*|S) = l/ 7 /)** 0) 

In equation (1), (5) represents the translational and rotational 
displacements and only includes those displacements which are 
associated with masses and rotary inertias. Other displacements 
are reduced out using static condensation. The spinning speed 
of the rotor is denoted by [G] is the gyroscopic matrix 
corresponding to the attached disks, (FI* represents the cou- 
pling forces on the rotor due to coupling to the housing while 
(F £ )* represents the external forces including the imbalance 
forces. Let 

(S|=1*]*I<7)a (2) 

in which ($1* is the modal matrix, normalized with respect to 
the mass matrix of the nonspinning rotor and |<7l* are the 
associated modal coordinates. 


_ Nomenclature 


['CO = structural damping ma- 
trix 

{Cl/ = coupling damping matrix 
[CM] = convergent matrix 
[D] = generalized damping ma- 
trix of rotor model 
e = imbalance of disk 
)£) = physical force 
[G] = gyroscopic matrix 
\H\ = physical coordinates of 
housing 

[V\] = identity matrix 
[K\ = structural stiffness matrix 
[XI/ = coupling stiffness matrix 
m = mass matrix 
( />) = generalized force (modal 
force) 

(<7) = generalized displacement 
(modal displacement) 


IQ] 

s 

(5) 

t 

T 

\y\ 


W|T» <*>d 

[ 0 ] 

M 

[ X A\] 


transition matrix 
physical displacement of 
rotor in bearing 
physical coordinates of 
rotor 
time 

time increment 
time at iT 

, generalized coordi- 
nates of rotor 
natural and damped nat- 
ural frequencies 
null matrix 
the rotor state matrix 
diagonal eigenvalue ma- 
trix 



<j> - spinning speed of rotor 
[4>] = normalized modal matrix 
f = damping ratio 


Brackets and Symbols 


[ 1 « 

v o = 

i i = 
i ]-' = 


Subscripts 

R = 
H = 
£ = 
I = 

x, r f z = 


square matrix 
diagonal matrix 
column matrix 
inverse of a matrix 
transpose of a matrix 


rotor 

housing 

external 

coupling 

X % T, Z direction 
/th time increment 
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Using the modal transformation (2), equation (1) is written 
in the form 

[q\*- (+1 W* I Q ) r + P ■ * f Q I * 

= (*l*(IP/)* + \Fe\r) ^ 

where 

rAJjr-lWJji- 

o>„ is the natural frequency of the nonspinning of the rotor. 
Adding a modal damping matrix [ X C\] to equation (3) yields 

+ PAxl*l?i*= [*)J(|P/)* + 1 F E I r > W 

where 

IZ>1 W = PCJ* - [*]J*[GW*J* (5) 

rco*-r2tto^i*. < 6 > 

(b) Housing. The equation of motion of the housing may 
be expressed as 

[M\ h \H) + \K\ h \H\ = [ F f ) H + \F e \ h (7) 

where [H\ represents the transverse displacement in the Y and 

Z-directions. In terms of the modal coordinates of the housing 
while uncoupled to the rotor, equation (7) takes the form (after 
adding a modal damping matrix) 

I^) w +[£»} w |<7) w +t x AJ//l9)w=[*]«((/ : '/)«+ \ f e ) h ) W 

where 

[L>] w =['C N ] w =[ x 2fu A 0 w (9) 

and 

( ,0 > 

(c) Coupling Forces. Forces due to linear coupling be- 
tween rotor and housing may include the direct and cross- 



turbines 


1! t Ml 



coupling forces due to seals, impeller, and turbine forces. The 
nonlinear coupling forces are taken as those due to the clear- 
ance between the rolling element bearing outer race and the 
housing. Figure 2 shows a model for the gap at each of the 
loosely supported bearings. 

The bearing force acting on the housing in the Y-direction 

is 


(Fq) y — 


Sy-Hy 


+ C G (S Y -H Y ) 

R<6 


(F C )y=C C {S Y -Hy) 

where the bearing support stiffness, K c > is shown in Fig. 3 and 
C c is the damping in bearing. Analogous equations hold for 
the Z-direction, and 


R = V<Sr-«i') 2 +(Sz-/*z) : 

where 

S y , S 2 are the physical displacements of rotor in the Y and 
Z directions, respectively, at the bearing location. 
H Yy H z are the physical displacements of housing in the Y 
and Z directions, respectively, at the bearing location. 

The bearing forces can further be written as 

(F C )y=K G (Sy-H Y )+C C (Sy-Hy) 

(Fc)z = A c (S z -//z) + Cc(5 z -///) 


where 


*g=*g(»-|) 
*c = 0 


Rs 5. 


(ID 

( 12 ) 

(13) 

(14) 


* Rfl. 1 Complex rotor system with flexible housing 


Method of Analysis 

A hybrid convolution method proposed by Noah et al.» 
(1986) and (1988), is utilized here for the analysis of the coupled 
rotor-housing system. The displacements of the housing are 
best represented by a convolution integral due to its accuracy 
advantage over the transition matrix formulation. Due to the 
presence of the unsymmetric gyroscopic terms, the equations 
of motion for the rotor are transformed to first order. The 
displacements of the rotor are expressed in terms of transition 
matrix for the motion of the rotor. 

The restriction of an impulse response convolution integral 
is that the equations of motion have to be decoupled. For a 
coupled subsystem, such as rotor with gyroscopic moments at 
the left-hand side of its equations, a transition matrix for- 
mulation has to be employed. However, the convolution in- 
tegral is more computationally efficient than transition matrix 
formulation. 

Transition Matrix for the Rotor. The equations of motion 
of the rotor, equation (4), are cast in first-order form, 

| w)* = [cr]*[u)*+ (15) 

where 


i z 




Fig. 2 Rotor and housing displacements at bearing location 
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If the spinning speed is a function of time, the gyroscopic 
term of equation (5) can be moved to the right-hand side of 
equation (4) to avoid solving a new eigenproblem at each time 
step. In this case, the coefficient matrices, given by equations 
(16) and (17), are replaced by 






_( (0] m \ 

M 


o 


[*]j£<MC]*[*],l<7l* y 
(17') 


(16') 


.) 


In this case, the equations of motion of the rotor becomes 
decoupled, and a convolution integral becomes a better selec- 
tion because it is more efficient than the transition matrix 
method. 

The solution of equation (15) can be written as 

(n(/)l* = e' lal ''(uol + j o e<'- r, '“>*|P(r))*dr (18) 


where 


Convolution Integral for Honsing. The decoupled form of 
the equations of motion of the housing allows expressing the 
nodal displacements of the housing by an integral. The integral 
formulation has the advantage of providing a closed-form 
expression as opposed to conventional numerical integration 
schemes. Also this representation allows dealing directly with 
diagonal modal matrices as opposed to coupled, first-order 
modal equations used with the Runge-Kutta method of the 
convolutional formulation matrices. This results in higher 
computational speed and accuracy of the convolution for- 
mulation. 

Based on equation (8), the housing generalized coordinates 
can be expressed as 
I q ( t) 1 H = [ V " fw " ; cos 0 [q{ 0 ))h 


+ 


+ 


r± 

Er= 


e f- "*sin 0 )^/ 


’\]({?(0))w+r«.i9(0)i«) 


-fwaU-r), 


sin u3 d (t 


-*] 


(P(r)Mr 


( 22 ) 


where ( q{ 0) ) , I ?(0) | are initial solutions, «*, are undamped 
and damped natural frequencies of housing, f is damping ratio 
of housing, (/>)„ = [«]£((f/ljr + \Fe\h). the generalized 
forces act on housing. 

The expression (22) is next written in discretized form, with 
the generalized forces taken as before as varying linearly within 
each time step, or 


(< 7 (' f+ i))w=r*~'"* , cos bf,. ,01 <7(0)1* 


[1 




sin ((<7(0)) w + a| <7(0)) w ) 


( sin bt ,+ , „„ cosfe/ l+ , } 

+ 031 

(ff(«..)lj#=r-e --, * I (« 005 bt i*i + b sin + 1 ) 1 <7(0) I w 


| sin bt i+[ + cos 


(( <7(0)1* 


+ a(d(0)|*) + 


-a sin bt^ , + b cos bt^ , 

btf + b 2 ) 


1 


fj . (19) 

*-o k ' 

Equation (19) is transition matrix and {u 0 l is the initial 
generalized coordinates. 

Equation (18) can be cast in a discretized form in which the 
generalized forces, ( P|*, are taken as varying linearly between 
time /, and / J + , so that 


iP(/> i = (P(/j) + t ~Y (|P(/i + i)i-|/ > (f<)l) 

where /,+ , = /, + Tand T is a small increment in time. The 
discretized form of equation (18) can then be shown to take 
the form 

\u(t,+ l )\ R = lQ(T)Uu(t i )\ ll 

+ (K2(T)1 - PADM*' (iW l* + Mi 1 


(/>('..!) l/i-tm 


/ ) 1 _ 


l*]iHlPUi+l)\*-lPVi)\K) 


( 20 ) 


where 

[Q(D] = e nal «. ( 21 ) 


a cos btj+i + b sin bt t + x 
btf + b 2 ) 




where 


a=fu n 
b = Ud 


EA M ~e- T EA, + A,^ 


EBn.,~e- T EB, + Bi+, 


EAi = e~ a,i ~'A i +e~'"'~ 1 Ai+ . . 

. +A 

EB, = e-" i -'B i +e-" t - 2 B i + . . 

. +B, 


(24) 

(25) 

(26) 

(27) 

(28) 

(29) 

(30) 


An.\ = Pi * i (a cos bt ,+ 1 + b sin bt it , ) 

_ (a 2 cos bt ,+ , + lab sin bt i+ , - b 2 cos bt ,+ , ) 

Ttf+b 2 ) 

Pjj^^Pj. e-'Tfa 2 cos btj + lab sin bt,-b ? cos bt,) 
TOt + b 2 ) 

~p l e-‘ T (a cos bt, + b sin bt,) 


(31) 
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t = P, + \( a sin bt,.i~b cos 

_ (a 2 sin bt,.]-2ab cos bt^^-b 1 %\n bt it \) 

7V+* 2 ) 

P,.\- p ±_ ~nT,g 2 sin bt-lab cos bt-trsm bt,) 

TUr + tr) 

-Pte~ ,T (a sin bt,- b cos bt,) (32) 

P is dement of vector (P| . 

The Coupling Forces. The coupling forces acting on the 
housing are given by 

[F,)h= l*l/(S) -[//)) + ICMIS ) - I »\ > (33) 

and those on the rotor are 

“ [F,) h 

where [tf], and Id, are the coupling stiffness and damping 
matrices, respectively. The coupling stiffness matrix includes 
the updated bilinear bearing stiffness at each iteration. 


(S3] = ^ ( n][*]M*l* ) 

,C, f -[T2 )\*] t h{K]A*)h - (72)[*]&[d/{*]« \ 
m = i -\mwSAC[A*\Hj 

and 

inj = [a]* 1 (aemi-r/o> ^y- -r/oj 

[ s sin bt^ i c _ 

[ btf + b 1 ) bia 1 - 


[72] 

173] = [ 


[ A O 


- a sin btj + x + b cos bt i+ 1 

b^ + b 2 ) 


acos 


b/ l4l + bsin br i4l 

b^ + b 2 ) J 


where 


The Computational Procedure. For the responses at time 
t = * f>l , the coupling forces between the housing and rotor 
are unknown. An iterative technique is used to calculate these 
forces. The following iterative loop is used in both the hybrid 


and Runge-Kutta methods: 

If the responses at / = tj are known and responses at t — f, + 1 
are desired, then: 

(a) Set |F;(/ I+1 )| = ( FiU,)]. 

(b) Calculate (P)* and {P)//. 

(c) Solve { w(/ f > i))* and \q(t^ i))//, l <7('< 4 1) )// by the hy- 

brid method (or the Runge-Kutta method). Calculate the Eu- 
clidean norm of (w(/ l+ 1))/? and i) )/*• 

(d) Transfer |u(/ J+ i)l* and {<?(/*+ i)lw. Iffa+iM## back 
to physical coordinates, then a new connection force (P/(f,> i) ) 
can be calculated. 

(e) If the Euclidean norms of two consecutive iterations 
of both rotor and housing are close enough then stop the 
iteration and go to step (0. otherwise go to (b). 

(0 Move a time increment and go to (a) until t reaches a 
specified time for terminating the run. 


The Convergence of Hybrid Method. Rewrite the coupling 
forces as the function of the generalized coordinates, equations 
(20), (23), and (24) can be arranged as the following form: 


- [cm * ( 


+ |c) (34) 


where k stands for *th iteration at a given time [CM] is 

the convergent matrix (Isaacson and Keller, 1966), |c) is vector 

with known values, and 


[CM] = 


where 

[SI] = [71] x 


/ [SI] [S2]\ 
y[S3] ]S4 ] J 

( 0 0 ^ 

l -[♦ jjtW / i *!* - [*] t r[cia*uJ 

( o 0 \ 

x y*]J[Al,[4>] w [*]1[C],[*]h J 


(35) 


AC=a cos b/, 4l + b sin bf,> t 


1 


n^+b 2 ) 


T^ + b 2 ) 

AS = a sin br 1+ , - b cos bt i+ 1 


(o 2 cos bti+ j + lab sin bf /4 j - b 2 cos b/,+ 1 ) 
(j 2 cos bf, + 2ab sin b^-^cos br,) 


r — — (j 2 sin bti+x-lab cos b^+i-b 2 sin bf J4 i) 

T^ + b 2 ) 

--•r 

+ — r- (a 1 sin bt:- lab cos bf ( -irsin b/,). 

n^+b 2 ) 


Equation (34) is in the form of the Jacobi iterative method 
(simultaneous iterations), sec (Isaacson and Keller, 1966). The 
convergent matrix, [CM] may not be Fixed at a given time 
/ = /, 4 , as in an ordinary linear system because matrix [AT] 7 may 
change during the iteration if the switch between two slope 
areas in Fig. 3 exists. [K], will be updated because the nonlinear 
bearing stiffness, K c . changes. In this case, there are two con- 
vergence rates during the iteration at that given time. 

Equation (34) will converge if and only if all eigenvalues of 
[ CM] are less than one in absolute value. An alternative way 
to check the convergence is that the iterative method will con- 
verge if, for any matrix norm (Euclidean, maximum, etc.), 

I [CM] I < 1. Uncommon cases occur where convergence cri- 
teria are met within one region in Fig. 3 and not the other. In 
that case, overall convergence at a given time j might occur 
depending on the pattern of movement from one region to the 
other. In any case, the convergence rate and size of the incre- 
ment Twill determine the outcome for acceptable accuracy. 

The rate of convergence, CR t is defined as the following: 

CR m - log p([CM]) 

where p {[CM]) is the spectral radius of [CM] defined by 
p([CA/])»maxlX,l 

where X, arc the eigenvalues of [CM]. 

Let the initial error be defined as the norm of the vector of 
the difference between the exact solution and an initial guess 
at the beginning of iteration. The number of iterations, y, 
required to reduce the initial error by the factor 10' m is in- 
versely proportional to CR and is defined as 

m 


Journal of Applied Mechanics 


R7 


SEPTEMBER 1990, Vol. 57 / 735 




Cc 


X 


ROTOR: Shaft ddsmeter OD = 7.26 x lO^m, CD = 0.0m 

Material: £ - 2.0684 x lO^Y/m* 

Joint length** 7.26 x 10 _1 m 
Rotor length^ 6.856 x 10" l m 

HOUSING: Houring/ rotor weight ratio * 6/1 

(£/)** 8.6947 x 
{EI)h= 17494 x \Q*Nm* 

4(0 = m, A',*^ 7.005 x 10* A/m C c * 0 0 
4(<) = K$y , - 8.7563 x 10*A/m 

f,{<) = m,e^coi^ Ksit~ 8 7 563 x 10 *A/m 
= m, e^ 1 rtn^i K jy> * 2.6269 x 10 t N/rn 
Fk(t) = m 1 e^ , ooa^t Ac * 8.7563 x 10 r A/m 
4(0 * TTViC+'tin+t S * 0.0127mm 

Rg. 4 The generic model 


Application and Discussion 

The present convolution method is applied to a modified 
version of a rotor model proposed by Davis et al. (1984). The 
model was proposed to represent a simplified generic mode 
of the SSME turbopumps. The parameters and coefficients of 
the present generic model (shown in Fig. 4) are given in Tables 
1 and 2. The imbalance forces are taken as shown on Fig. 4. 

Response of the generic rotor-housing model at the bearing 
location is determined for the hypothetical startup-shutdown 
case shown in Fig. 5. For various time increments, A/, a com- 
parison is made of the errors and computer CPU time (on a 
VAX 8650) between the results obtained using Runge-Kutta 
fourth-order method and those using the present convolution 

method. , , . 

Comparisons of the CPU time, are made with comparable 
errors of the results of using both methods. The errors is 
measured at time t = 0.01 seconds and defined as the following 


r I v„-v\ 

Error = Tv 

where V n is the exact solution of the rotor and housing dis- 
placements in vector form and V represents the displacements 
of rotor and housing at the corresponding time in vector form. 

The “exact” solution is calculated using the hybrid con- 
volution method with a very small time increment, 2 x 10 
seconds and tolerance = 1 x 10* i: . The tolerance is the 
absolute value of the difference of the displacement and ve- 
locity norms of rotor and housing described in step (e) of the 
computational procedure. 

As shown in Figs. 6 and 7, with the tolerance equals to 1 
x 10“ 8 seconds, the hybrid convolution method is faster and 
more accurate. For more meaningful comparison, the accur- 
acies of the Runge-Kutta and convolution methods are made 
closer by reducing the allowable tolerance for the Runge-Kutta 
to 1 x 10” 12 seconds and 1 x 10“ 14 seconds. The CPU time 
for the Runge-Kutta increases quickly from 1.42 times of the 
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Table i PimMltri ol the generic rood#* 


Disk 

Mu. 1 

Imbalance 

Moment of inertia ( A f m J ) 

{*9) 

(mm) 

Polar 

Diametrical 1 

r— l 

i 

lL 

4.5358 

6.8124 

9.0175 

0.0908 

0.0906 

0.0908 

3.1635 x 10'* 
2.1806 x 10" 3 
1.1363 x lO' 1 

1.7064 x 10* s 

1.1411 x 10‘ 3 

5.7373 x lO' 1 
— 


Table 2 CoaHtctonta of the generic model 


Disk 

Kyy-i Arr 
( N/m) 

A'xv, Ay* 
(N/m) 

Cyy i Ciz 
(JV $tcjm) 

1 

1.7513 x 10 T 

2.6269 x 10* 

1.1033 x 10 4 

2 

-1.5586 x 10* 

-7.0050 x 10" 

1.7513 x 10 1 

3 

1.5236 x 10* 

3.3274 x 10* 

4.0279 x 10* 



Side force {N) 

Disk 

Y-direction 

Z-direction 

1 

2.2019 x 10"* 4* 

-1.2482 x UT 4 4> 3 

2 

2.7210 x lO -4 ^ 

3.2387 x 10‘ 4 & 

3 

1.0582 x 10 " 4 ? 

0 


(rad/aec) 



hybrid method’s CPU to 4.23 times (with a time increment of 
2 x 10“ 5 ). For the closest accuracies, the CPU time of the 
Runge-Kutta is 4.23 times that of the hybrid method. The 
hybrid convolution method is also more robust than the Runge- 
Kutta method. The Runge-Kutta algorithm failed to converge 
for time increments greater than 2 x 10 5 seconds. However, 
the hybrid method will diverge when the increment used is 
larger than 3.3 x 10“ 5 seconds. 

Concluding Remarks 

The convolution integral and transition matrix methods are 
represented in closed form. They generate less roundoff errors 
and require less numerical computation time in comparison 
with the Runge-Kutta fourth-order method. Moreover, they 
are more robust than the Runge-Kutta method since they will 
converge for larger time-step size. 

The hybrid convolution approach developed in this study is 
shown to provide an efficient and accurate closed-form integral 
formulation for determining the transient response of linear 
systems coupled through local nonlinearities. A typical appli- 
cation in which the present method proved quite effective is 
the determination of the transient response of a generic model 
of the high pressure oxygen turbopump (HPOTP) of a space 
shuttle main engine (SSME) in presence of bearing clearances, 
constituting the local nonlinearities. Substantial savings in 
computation time were achieved as compared with direct nu- 
merical integration techniques. 

The use of the transition matrix allows the representation 
of rotors involving skew-symmetric matrices of gyroscopic 
loads or other nonconservative systems with general velocity 
coefficient matrices. A convolution integral would represent 
quite effectively other systems with classical modes, such as 
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0 T -* 

LEGEND 

<*•*- hybrid (T0c.-1.E-8) 

_+■ RUNGE KUTTA (T0L.-1.E-B) 

05- - *- RUNGE KUTTA (T0L.-1 -E- 1 2) 

-O- - RUNGE KUTTA (T0L.-1.E-M) + 



TIME STEP SIZE (X 10000 SECOND) 

FJy. 6 Accuracy of hybrid mathod and Runo^Kutla mathoda 



TIME STEP SIZE (X 10000 SECOND) 


Ftg. 7 CPU tima ol hybrid mathod and Runga-Kutta mathoda 


the housing of the HPOTP or other nonrouting, proportion- 
ally damped structures. The convolution formulation allows 
accommodating with ease changes in the nonlinear or linear 
coupling parameters among the various linear subsystems in- 
volved. 

Possible improvement of the method could be achieved 
through other alternative or optimization for the iteration pro- 
cedure utilized in this study. The methods described in this 
study is believed to be capable of handling any type of local 
nonlinearity. For an impact case at the bearing, for example, 
iteration over the contact force could be replaced by iteration 
to esublish contact/no-contact at each increment and then 


reverse the impact velocity. Applications of the present method 
to systems with various types of local nonlinearities could be 
worthwhile. 
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